- Email: [email protected]
A Scientific and Statistical Analysis of Accelerated Aging for Pharmaceuticals. Part 1: Accuracy of Fitting Methods
RESEARCH ARTICLE – Drug Discovery-Development Interface
A Scientific and Statistical Analysis of Accelerated Aging for Pharmaceuticals. Part 1: Accuracy of Fitting Methods KENNETH C. WATERMAN, JON T. SWANSON, BLAKE L. LIPPOLD FreeThink Technologies, Inc., East Lyme, Connecticut Received 18 February 2014; revised 31 May 2014; accepted 4 June 2014 Published online 8 July 2014 in Wiley Online Library (wileyonlinelibrary.com). DOI 10.1002/jps.24075 ABSTRACT: Three competing mathematical fitting models (a point-by-point estimation method, a linear fit method, and an isoconversion method) of chemical stability (related substance growth) when using high temperature data to predict room temperature shelf-life were employed in a detailed comparison. In each case, complex degradant formation behavior was analyzed by both exponential and linear forms of the Arrhenius equation. A hypothetical reaction was used where a drug (A) degrades to a primary degradant (B), which in turn degrades to a secondary degradation product (C). Calculated data with the fitting models were compared with the projected roomtemperature shelf-lives of B and C, using one to four time points (in addition to the origin) for each of three accelerated temperatures. Isoconversion methods were found to provide more accurate estimates of shelf-life at ambient conditions. Of the methods for estimating isoconversion, bracketing the specification limit at each condition produced the best estimates and was considerably more accurate than when extrapolation was required. Good estimates of isoconversion produced similar shelf-life estimates fitting either linear or nonlinear forms of the Arrhenius equation, whereas poor isoconversion estimates favored one method or the other depending on which condition C 2014 Wiley Periodicals, Inc. and the American Pharmacists Association J Pharm Sci 103:3000–3006, 2014 was most in error. Keywords: isoconversion; accelerated stability; shelf-life; Arrhenius fitting; kinetics; stability; chemical stability; mathematical model; in silico modeling
INTRODUCTION Determining the stability of drug substances and drug products is one of the critical, and at times, gating elements in the drug development process. In recent years, new protocols and mathematical models have been applied to the use of accelerated processes for rapidly determining the chemical stability of APIs and drug products both in solution and solid state.1–10 With the increased use of these models, it is appropriate to consider the implications of different methodologies in terms of both the statistics and science involved. In Part 1 of this detailed examination, the ability of different models to accurately estimate shelf-life is determined on a theoretical basis. Part 2 is focused on different statistical modeling methods for estimating the precision of projections.11 Finally, Part 3 examines the calculation options when there is little or no measured degradation.12 While the current paper examines what mathematical model best fits data without error bars, the second paper builds from that model to include statistical error propagation. The third paper evaluates the implications of different options within the model and error bar propagation of the first two papers for cases where there is little degradation occurring with respect to the error bars for the measurements. While previous examinations of accelerated stability modeling were based largely on loss of API potency with time,13–19 currently, increase in related substance (degradant) growth is more commonly the limiting factor on setting shelf-lives. This has implications in the modeling process for three reasons. First, the specification limit for the loss of potency (usually between 5% and 10% loss) is typically higher than the specification limit for Correspondence to: Kenneth C. Waterman (Telephone: +860-287-4253; E-mail: [email protected]) Journal of Pharmaceutical Sciences, Vol. 103, 3000–3006 (2014)
C 2014 Wiley Periodicals, Inc. and the American Pharmacists Association
3000
increase in the level of a degradant (typically 0.2%–1%). Second, the precision of the measurements for related substance growth is generally much tighter than that for potency loss. Third, in secondary degradation processes, potency loss can show simple kinetics even when individual degradant formation is complex. In classic chemical transformations, processes can be described by a reaction order related to the species involved in the rate-limiting step of the degradation process. An assumption is often made that the kinetics are pseudo-zero order because the extent of degradation is typically low (0.2%–1.0%) such that a straight line can be fitted to measured degradation versus time points. The rate under the measurement conditions will be the slope of that line; that is, the initial rate of degradant formation is often assumed to determine the rate constant. In solid-state and complex solution-state pharmaceutical systems, there are a number of reasons that this simple reaction order (and initial rate based) model does not describe the situation adequately, especially when focused on degradant formation rather than loss of API. First, degradants can undergo further degradation in secondary processes. This affects the primary degradant (rate slows as reaction proceeds) and secondary degradant (a lag followed by growth) formation kinetics. Second, diffusion processes can be rate-limiting steps, especially in solid-state chemical reactions. The kinetics of such reactions does not follow the simple reaction order model. Third, degradation products can catalyze further degradation in an autocatalytic process. This results in a lag followed by more rapid growth kinetics. A similar kinetic trend is seen in cases where an inhibitor (e.g., antioxidant) is consumed as the reaction proceeds. Finally, many solids are heterogeneous such that different nonequilibrating forms of the API react at different rates. Examples of this are drug molecules at interfaces versus in the bulk, and molecules at defect sites versus in crystalline domains. In our examination of calculation options, one of the evaluation factors for different
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
RESEARCH ARTICLE – Drug Discovery-Development Interface
methods will be the ability to accurately estimate behavior with complex kinetics. In experience with many solid drug products, we have seen that more than half do not follow a simple kinetic model. One of the processes for stability assessments evaluated here that has seen recent increased use in pharmaceutical applications is the isoconversion method.2–4 With isoconversion, rather than a focus on rate constants, the process involves an emphasis on the time hit the specification limit (the isoconversion time). For processes that can be characterized by a rate constant, the isoconversion time will be the specification limit divided by that rate constant. As a simple example for complex kinetics, we consider a process where a drug (A) goes to a primary degradation product (B), which undergoes further degradation to C, with the rate constants for the two processes corresponding to k1 and k2 : A→B→C
(1)
The amount of B and C as a function of time will correspond to Eqs. 2 and 3 (derived from Ref. 20): Bt = A0 Ct = A0 1 +
−k t k1 e 1 − e−k2 t k2 − k1
(2)
−k t 1 k2 e 1 − k1 e−k2 t k1 − k2
(3)
The behavior of B and C as a function of time is shown in Figure 1. As seen in the figure, B will rise with a decreasing rate, whereas C rises after a lag (as B builds up). Since the equations governing these species are known, one can compare the results from an exact solution of the problem to those obtained using the different methods for using accelerated data to estimate ambient stability. The kinetic forms reflected in Eqs. 2 and 3
3001
will be representative of many of the types of observed non-zero order behaviors in drug products. The dependence of solid-state chemical reactions on temperature and RH has been shown to follow a modified Arrhenius equation which can be expressed in linear (Eq. (4) and exponential (Eq. (5) forms. ln k = ln A − Ea /(RT) + B(RH)
(4)
k = A e(−Ea /RT)+B(RH)
(5)
where k is a rate for degradant formation, T is the temperature (K), and RH is the equilibrium relative humidity to which the sample is exposed (as opposed to the external RH for a packaged product).1,2 Varying T and RH to determine the rate at each condition, then determining the best fit of the variables A, Ea , and B allows a mathematical model to be built which can be used to calculate the rate (and corresponding shelf-life) at extrapolated or interpolated T and RH conditions. In real applications of accelerated aging, the RH dependence complicates predictions of ambient behavior yet cannot be ignored. Because the form of the equations remains the same with or without the humidity term, the calculations herein assume a constant RH throughout with the implication that the preferred methodologies for accurate predictions of shelf-lives will be the same when the RH term is added. The added complexity of the modified Arrhenius equations, however, will exacerbate any inaccuracies intrinsic in a method. For purposes of illustration, we are using a system where the true kinetic form of the degradant formation is known. In practice, this is often not the case. One method for determining shelf-life when the kinetic form is not known involves calculation of the isoconversion time; that is, the time when the degradant reaches its specification limit at each condition. At the specification limit, the isoconversion time will be proportional to the (isoconversion) rate. This allows the isoconversion time to be used to calculate the Arrhenius parameters at the specification limit, without knowing the form of the rate equation. The shelf-life at ambient conditions is the isoconversion time at those corresponding conditions. One can illustrate the relationship mathematically for several simple cases. In measuring degradant formation (d) for a zero-order or pseudo-zero order reaction with a rate constant k, degradant formation will be proportional to time (t) as shown in Eq. (6): d = −kt
(6)
The time at which the degradant reaches the specification limit is defined as the isoconversion time (iso). Since d at each isoconversion time is the same at all temperatures (it is the specification limit), the following will be true: kT1 iso1 = kT2 iso2 Figure 1. Model drug behavior used in the present method evaluation showing the 60◦ C curves. Drug substance A degrades to a primary degradant B (blue curve), which in turn degrades to a secondary degradant C (orange curve). The graph shows B and C as a function of time and typical specification limits of 0.2% and 0.5% for degradant formation. DOI 10.1002/jps.24075
(7)
Rearranging this equation gives the result that the ratio of isoconversion times is inversely proportional to the rate constants: kT iso2 = 1 iso1 kT2
(8)
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
3002
RESEARCH ARTICLE – Drug Discovery-Development Interface
One can substitute into the rate equation (Eq. (6) the specification limit, the rate constants (as determined by the Arrhenius equation) and isoconversion times at two temperatures to determine the activation energy (Eqs. 9–12). spec = A e−Ea /RT1 iso1
(9)
spec = A e−Ea /RT2 iso2
(10)
E iso2 − a =e R iso1
1 T1
R ln Ea = −
1 T1
− T1
iso2 iso1
−
2
1 T2
(11)
(12)
Once the activation energy is determined, one can calculate the ratio of rates between a known temperature and any other temperature. Determining the isoconversion times at two temperatures is therefore the minimum requirement to determine the shelf-life at ambient conditions. The true rate does not need to be known. This can be extended to more complex kinetics. For example, first order kinetics follows Eq. (13): d = [A]0 1 − e−kt
(13)
d − [A]0 = −e−kt [A]0
(14)
or
Taking the natural logarithm gives Eq. (15). − ln
d − [A]0 [A]0
= −kt
(15)
This is the same as Eq.
(6) with d replaced by the expresd−[A]0 0 , and − ln will be the same at differsion − ln d−[A] [A] [A] 0
0
ent temperatures at the specification limit. Therefore, Eq. (12) again allows one to relate the isoconversion time at an accelerated temperature to the isoconversion time at an ambient condition (shelf-life). In fact, no matter what the form of the degradant formation kinetics, using isoconversion will maintain this relationship, assuming the reaction mechanism to form the degradant does not change with temperature. This means that the shelf-life can be extrapolated from accelerated conditions without knowing the true kinetic form for degradant formation as long as the isoconversion times are known at two or more temperatures. For the present evaluation, rather than illustrate a more complex example using equations, several methods of shelf-life prediction will be compared with the value one would obtain if the true kinetic form was known. Three methods for estimating degradant formation rates are compared with respect to their ability to provide accurate estimations of shelf-life from accelerated data.
First, the amount of degradant at each time point is considered separately with respect to the initial amount of degradant. Effectively, each measured value at time ti and temperature Ti will result in a rate constant ki based on the change in the amount of degradant divided by time. These values can be used to determine a set of independent rates based on all time points and conditions. Second, each temperature condition is used to determine a best fit, zero-order rate constant for that condition. This is effectively a rate based on the initial slope for degradant formation. Multiple time points contribute to the rate estimation in each condition. Third, an isoconversion method is employed. In practice, the isoconversion time point will not be known a priori, and multiple measurements of degradant formation as a function of time will be made. These measurements can be used to estimate the isoconversion time at each condition. How best to use these data for the estimation of isoconversion will be examined using three methods. (1) A line will be drawn through the data to find the intercept with the specification limit. (2) The two closest points to the specification limit will be used to either interpolate or extrapolate to the specification limit. (3) The single point closest to the specification limit will be used to define a rate using k = d/t. Each of the methods results in a different rate or set of rates at a given condition, which will have implications when these rates are used with different temperatures to predict shelflife, particularly when the extent of reaction is different at the different temperatures. In those cases where ideal linear (zero order) behavior, without noise, is observed (or when the rate can be characterized effectively based on the initial slope), all the methods will result in the same rate constant. In our examination of drug product stability, over half the cases studied (n > 200) do not show such ideal behavior. In addition, since most pharmaceutical stability studies used in product development involve very few time points, it can be difficult to validate linearity in a given system. We will examine cases where there is apparent linearity, but where the degradant level remains well below the specification limit to compare how accurately the various treatments predict ambient behavior in such cases. Two methods are examined for determining the temperature dependence of degradant formation and thereby estimating the shelf-life at ambient conditions. Each assumes the temperature dependence can be modeled by the Arrhenius equation. (1) The equation can be transformed to a linear form (Eq. (4) and parameters determined by a linear least squares method. (2) The equation (Eq. (5) can be used directly and parameters determined using a nonlinear least-squares method. In either case, one is minimizing an objective function, typically the sum of the squared error between the actual and predicted values. If the values fit the theoretical model perfectly, both methods will give the same results. With real data, however, this situation will not be true, and the two methods will produce different results because of the difference in the objective function. When fitting Arrhenius data, the relative leverage of the higher temperature points will be greater for the nonlinear fit, since the predicted value depends exponentially on the temperature. Loss of potency can be particularly difficult in some cases for Arrhenius fitting since multiple pathways for loss of API can each have different activation parameters; however, this issue is not evaluated in this study. Similarly, differences in the
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
DOI 10.1002/jps.24075
3003
RESEARCH ARTICLE – Drug Discovery-Development Interface
Table 1.
Primary Degradant-Based Shelf-Life with Different Modelsa Shelf-Life (yrs) 25◦ C
Method Specification Limit 0.2%
Exactb Single constant time point (14 days)c Individual time pointsd Line through all pointse Single point at isoconversionf Linear, all-point isoconversiong Bracketed isoconversionh Closest point isoconversioni
Specification Limit 0.5%
Linear
Nonlinear
Linear
NonLinear
1.43 0.56 0.62 0.29 1.43 12.35 1.36 1.23
1.43 0.40 0.54 0.19 1.43 467.94 1.26 1.25
4.45 1.40 1.56 0.71 4.45 1.40 3.19 2.76
4.45 1.00 1.36 0.49 4.44 1.16 4.39 2.63
a Shelf-life at 25◦ C estimated from accelerated conditions of 50◦ C, 60◦ C, and 70◦ C for a hypothetical drug product with the shelf-life being limited by formation of a primary degradant that can degrade to a secondary degradation product (k1 and k2 at 50◦ C equal 0.000113 and 0.011250%/d). Time points at each temperature, in addition to the origin, were 3, 7, 14, and 28 days b Refers to the true value based on the rate constants. c Based on linear rate constants using only one time point at each temperature. d Based on four linear rate constants at each temperature. e Based on using linear rate constants through four time points at each temperature. f Using the exact time points corresponding to the intersection of the degradant formation curves with the specification limit. g Estimating isoconversion times by linear fitting of all four time points at each condition to determine the intersection with the specification limit. h Estimating isoconversion times using the two time points closest to the specification limit at each temperature (or extrapolating from the last two time points when the specification limit is greater than the values at that temperature). i Estimating isoconversion times using the point (one of four) closest to the specification limit at each temperature.
activation parameters for the two rate constants in our A to B to C scenario are not considered here. In the present investigation, the focus is on the impact of these different calculation methods on the accuracy of predictions using perfect (i.e., no error) data, but with complex (but common) kinetics. It can be asserted that the ability of a method to provide accurate predictions for these types of kinetics will extend to other complex (or simple) kinetic schemes. We assume that the true kinetic form is unknown in the investigation, which is generally the case in practice, since it is rare that stability studies are carried out with enough time points and duration to truly determine the kinetic form. Even when kinetics appears linear over a range, curvature at greater conversions (nearer to the specification limit) may significantly affect the estimations.
METHODS Many of the calculations were carried out using the R statistical package.21 The R statistical package contains a number of subroutines used herein, including the “lm function,” which is a linear fitting module, and “nls function,” which is a nonlinear, least squares fitting routine. Both these routines are available as part of this open-source package. Values of the percent degradant as a function of time were calculated using Eqs. 2 and 3. The rate constants used met the following three criteria: the isoconversion time with an 0.2% specification limit and a temperature of 50◦ C is within the 28 day window of the hypothetical accelerated stability study; the isoconversion time with an 0.5% specification limit and a temperature of 50◦ C exceeds that window (approximately 60 days); the shelf-life at 25◦ C with a 0.5% specification limit is approximately 4 years (activation energy of 25 kcal/mol). Rates at 60◦ C and 70◦ C were determined assuming the 25 kcal/mol activation energy. For the calculations of B in the A to B to C scheme, k1 and k2 at 50◦ C equal 0.000113 and 0.011250%/d, respectively; for calculations of C, k1 , and k2 at 50◦ C equal 0.000112 and 0.090000%/d, DOI 10.1002/jps.24075
respectively. The points used for the calculations at each temperature were 0, 3, 7, 14, and 28 days. Three methods were used to calculate the rate constants used for determining the Arrhenius parameters. The first method calculates the rate of each individual time point as k = d/t. This method was used for calculations with a single constant time point (d was the value determined at a time of 14 days), a single measured point closest to isoconversion, and multiple individual time points. The second method fit the points to a line using the lm function of the R statistical package.21 The third method calculated k as k = specification limit/isoconversion time. This method was used for the single point at isoconversion, the linear fit to isoconversion, and the bracketed isoconversion. In the latter two cases, the isoconversion time was first calculated by determining the intercept of the best fit line with the line d = specification limit. The lm function was used to fit the points (all or the bracketing points). The Arrhenius parameters were determined by fitting Eq. (16) (linear fit) using the lm function or Eq. (17) (nonlinear fit) using the nls function of the R statistical software package. ln k = ln A −
Ea RT
k = Ae−Ea /RT
(16)
(17)
The Arrhenius parameters were used to calculate k at 25◦ C, which was then used to calculate the shelf-life as specification limit/k.
RESULTS AND DISCUSSION Table 1 summarizes the results of the calculations using different methods for determining the extrapolated shelf-life of the B degradant in an A to B to C scheme. Results reflect the use of Eq. (2) to calculate degradant formation. The shelf-life was
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
3004
RESEARCH ARTICLE – Drug Discovery-Development Interface
Table 2.
Secondary Degradant-Based Shelf-Life with Different Modelsa Shelf-life (yrs) 25◦ C
Method Specification Limit 0.2%
Exactb Single constant time point (14 days)c Individual time pointsd Line through all pointse Single point at isoconversionf Linear, all-point isoconversiong Bracketed isoconversionh Closest point isoconversioni
Specification Limit 0.5%
Linear
Nonlinear
Linear
Nonlinear
2.02 8.29 16.64 3.29 2.02 2.75 2.06 2.14
2.02 4.24 6.84 2.15 2.02 2.05 2.05 3.12
4.01 20.74 41.61 8.21 4.01 7.56 4.78 7.41
4.01 10.60 17.11 5.36 4.01 5.23 4.25 6.56
a Shelf-life at 25◦ C estimated from accelerated conditions of 50◦ C, 60◦ C, and 70◦ C for a hypothetical drug product with the shelf-life being limited by the formation of a secondary degradant (k1 and k2 at 50◦ C equal 0.000112 and 0.090000%/d). Time points at each temperature, in addition to the origin, were 3, 7, 14, and 28 days b Refers to the true value based on the rate constants. c Based on linear rate constants using only one time point at each temperature. d Based on four linear rate constants at each temperature. e Based on using linear rate constants through four time points at each temperature. f Using the exact time points corresponding to the intersection of the degradant formation curves with the specification limit. g Estimating isoconversion times by linear fitting of all four time points at each condition to determine the intersection with the specification limit. h Estimating isoconversion times using the two time points closest to the specification limit at each temperature (or extrapolating from the last two time points when the specification limit is greater than the values at that temperature). i Estimating isoconversion times using the point (one of four) closest to the specification limit at each temperature.
then calculated for two different specification limits, with different modeling methods using both the linear and nonlinear versions of the Arrhenius equation. Similarly, Table 2 provides the results for the C degradant. Formation of a primary degradant that undergoes further degradation to a secondary degradant is a common example of complex degradation kinetics. Traditionally, stability projections using accelerated data have assumed linearity (zero order or pseudo zero order kinetics, or use of initial slope). As seen in the results for the primary degradant (Table 1), whether using a single time point or multiple points, this assumption leads to significantly inaccurate estimations of room temperature behavior. With formation of a primary degradant, the linearity assumption underestimates the ambient shelf-life. This occurs because the lower temperature conditions do not show the downward curvature to as great an extent as with the higher temperatures. This leads to a lower apparent temperature slope (activation energy) and a rate higher than the true rate at ambient conditions. With formation of a secondary degradant, there is a lag before growth. The results with an assumption of linearity (Table 2) are again quite inaccurate. In this case, the lag phase is not yet exhausted at lower temperatures. Since higher temperatures show greater reactivity, the linear approximation results in an apparent higher activation energy and an overestimation of the shelf-life at ambient conditions. While the primary degradant would give an estimation that is overly conservative, with a secondary degradant, the linear approximation could result in product failures on storage. As long as the degradant formation (whether primary or secondary) is assumed to be linear (zero order or pseudo zero order), the estimations to ambient conditions from high temperatures are not in general improved by using a nonlinear least squares fitting to the Arrhenius equation. The rates themselves at the accelerated conditions do not accurately represent the behavior to the specification limit. The form of the Arrhenius equation that is fit cannot fix this problem.
The situation is significantly improved when using isoconversion times. With the exact isoconversion times used for each condition, the projections to ambient stability are exactly correct for both primary and secondary degradants (Tables 1 and 2) as anticipated in the introduction. The practical challenge for using accelerated aging in such systems is how to determine the isoconversion times when they are not known a priori. Three different approaches at two specification limits were evaluated in this study. There are two reasons for examining multiple specification limits. Whereas with the zero-order assumption one has the same rate at all times, this is not true for complex kinetics. Because the specification limit is the time at which the degradant reaches the specification limit, the effective rate will depend on the specification limit and therefore the Arrhenius parameters (A, Ea ) will also depend on the specification limit. The two cases also illustrate another important difference. At the lower specification limit, all isoconversion times can be interpolated. At the higher specification limit, extrapolation is needed for the 50◦ C point. This allows comparison of the accuracy of interpolated and extrapolated data. The most obvious method for estimating isoconversion is to use a linear fit of the time versus percent degradant data to determine the time to reach the specification limit. The isoconversion rate would then be the difference between the specification limit and the initial degradant level divided by the isoconversion time. This is effectively the same as the linear kinetics assumption (zero or pseudo zero order), except that here the rate propagated to the Arrhenius fitting is the rate to the specification limit. This method creates a significant error when the fitted line is a poor estimation of the true isoconversion. This will happen when most of the points used to fit the line are well beyond isoconversion or extrapolation is required to reach the isoconversion time. Whether the isoconversion time is overestimated or underestimated depends on the curvature of the degradant time plot and whether the isoconversion time is being interpolated or extrapolated. These factors also influence
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
DOI 10.1002/jps.24075
RESEARCH ARTICLE – Drug Discovery-Development Interface
Figure 2. When including points well past isoconversion (illustrated here with degradant B at 70◦ C), a linear fit to the data may give a nonzero intercept. In this example, this results in underestimating the isoconversion time (estimated isoconversion time of 0.56 days shown in red versus the true isoconversion time of 2.04 days, shown in green).
whether it is the higher or lower temperature points which have the least accurate estimate of the isoconversion time. A second strategy to estimate isoconversion evaluated here is to use only the measured points that bracket the specification limit (or extend the closest two for extrapolations), rather than use all the points measured. This better accounts for the fact that the rate varies with time. As can be seen in Tables 1 and 2, the estimations with this process are significantly better for both the primary and secondary degradants. While not perfect, the “bracketed” approach provides a reasonable estimation of the shelf-lives in all cases. Comparing the predictions from the different methods to estimate isoconversion in the examples illustrates several key points. First, using all four time points works particularly poorly when estimating the shelf-life of the primary degradant at a 0.2% specification limit. This is the result of a single poor estimation—the isoconversion time at 70◦ C (Fig. 2). When points well past the isoconversion time are used, as in this ex-
3005
ample, the isoconversion time is underestimated and a serious overestimation of the shelf-life results. This problem is easily remedied by not including points significantly past the isoconversion time. A more difficult problem occurs when the study is not carried to isoconversion, as happens at 50◦ C with the 0.5% specification limit (Fig. 3). In this example, even though the points measured show nearly perfect linearity (R2 = 0.998), the nonlinearity of the true kinetic form is observed by the time isoconversion is reached. Linearity observed over the range of the experiment is insufficient to assume zero-order kinetics. Degradant formation must remain linear out to the isoconversion time. The result is a significant underestimation of the isoconversion time at 50◦ C, which results in a corresponding underestimation of the shelf-life. In general, when kinetics is complex, extrapolated values can be much less accurate than interpolated values, which are reflected in the estimate of the shelf-life. Bracketing the isoconversion time works well in these examples because the data are exact. In real studies, one must balance minimizing the error that results with complex kinetics when a linear estimation is used, with the error inherent in measurement of the data points. In the former case, focusing on the isoconversion time will reduce the error. In the latter case, more measurements will provide a more accurate estimate of the true degradant levels as a function of time. The commercial software ASAPprime incorporates an algorithm to look at how many points are appropriate to use for the estimation in the event the overall fit cannot be approximated as linear. The final estimation method for isoconversion uses just the point closest to the specification limit to define a rate. As seen in Tables 1 and 2, the estimations from this method are not as good as with the “bracketed” process; however, they are better than most other methods. The fact that a single point near isoconversion provides a better estimation than most of the other methods illustrates the importance of accurately estimating the isoconversion time for complex kinetics. This is not a method we would recommend, as per the preceding discussion, it benefits in this simple example from the fact the data points exactly fit the kinetic scheme. R
Figure 3. When measurements are not carried out to isoconversion (illustrated here with degradant B at 50◦ C), curvature may not be apparent (R2 for the fitted line is 0.998). In this example, the isoconversion determined by a linear fit (estimated isoconversion time of 51.6 days shown in red versus the true isoconversion time of 62.0 days, shown in green), emphasizing the desirability of carrying out stability studies to the specification limit at each condition. DOI 10.1002/jps.24075
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
3006
RESEARCH ARTICLE – Drug Discovery-Development Interface
When nonlinear fitting to the Arrhenius equation is employed with the isoconversion process, the results in most cases are not significantly different from those with the linear fitting. When the estimation method is exactly correct, the linear and nonlinear fits produce the same result. When the estimation method is not exact, since the nonlinear fit leverages the higher temperature points more heavily than the linear fit, it will tend to give better estimations when the high temperature estimations of isoconversion are closer to the true isoconversion values. The opposite will be the case when the low temperature values are better estimations. In the examples used, the higher temperature conditions required only interpolation to determine the isoconversion times. The lower temperature conditions in some cases required extrapolation. The result is that in some cases the nonlinear fit appears to work better. The more significant difference between these methods relates to the estimated uncertainties of prediction produced by these methods, which is outside the scope of the present paper (discussed in the subsequent papers11,12 ). It is possible to include a weighting factor onto either linear or nonlinear least-squares fitting methods. In particular, since nonlinear (exponential) fitting biases toward the high temperature points, it should be possible to counter this by weighting the lower temperature points more heavily. Ultimately, this opposes the argument that using the “true” form of the Arrhenius equation intrinsically should be more accurate.14 While it is likely that with appropriate weighting nonlinear methods could approximate the linear methods; however, since the latter are computationally less complex, it makes more sense to simply use linear fitting. Moreover, any weighting will also play a part in the statistics for the uncertainty of the prediction, which is discussed in the subsequent paper in this series.11
CONCLUSIONS Traditional (i.e., linear rate constants estimated by initial slope) and isoconversion methods were compared for the accuracy of predicted shelf-lives from accelerated data when the kinetics of primary and secondary degradant formation was complex. The results show clearly that traditional stability assumptions provide poor estimates. This is true whether a single or multiple time points are used, and whether a linear or nonlinear Arrhenius fitting is used. With isoconversion calculations, shelf-life predictions are directly related to how good an estimate of isoconversion is achieved at each condition. Of the methods examined, bracketing of the specification limit (extending the last two time points for extrapolations) provides the best predictions. The best results are obtained when the isoconversion estimates are obtained by interpolation. Even relatively short extrapolations can introduce considerable error. In this study, linear and nonlinear fitting (parameter estimation) of the Arrhenius equation provide similar accuracy for ambient shelf-life predictions. A key point is that if the isoconversion times are accurate, the two approaches produce the same result. It is not possible to decide whether linear or nonlinear fitting provides better accuracy since the results depend on which condition is in error. A proper discussion of these two approaches requires also considering the estimate of uncertainty they produce (as discussed in the second paper in this series11 ). Conflict of interest: The authors have a commercial stake in the software ASAPprime , which broadly uses accelerated aging data with some of the methods described in the manuscript. R
REFERENCES 1. Waterman KC, Adami RC. 2005. Accelerated aging: Prediction of chemical stability of pharmaceuticals. Inter J Pharm 293:101–125. 2. Waterman KC, Carella AJ, Gumkowski MJ, Lukulay P, MacDonald BC, Roy MC, Shamblin SL. 2007. Improved protocol and data analysis for accelerated shelf-life estimation of solid dosage forms. Pharm Res 24(4):780–790. 3. Waterman KC, Colgan ST. 2008. A science-based approach to setting expiry dating for solid drug products. Regul Rapporteur 5(7):9–14. 4. Waterman KC. 2009. Understanding and predicting pharmaceutical product shelf-life. In: Huynh-Ba K; Ed. Handbook of stability testing in pharmaceutical development. New York: Springer, pp 115– 135. 5. Waterman KC, MacDonald BC. 2010. Package selection for moisture protection for solid, oral drug products. J Pharm Sci 99:4437–4452. 6. Waterman KC. 2011. The application of the accelerated stability assessment program (ASAP) to quality by design (QbD) for drug product stability. AAPS PharmSciTech 12:932–937. 7. Colgan ST, Watson TJ, Whipple RD, Nosal R, Beaman JV, De Antonis DM. 2012. The application of science- and risk-based concepts to drug substance stability strategies. J Pharm Innov 7(3–4):205–213. 8. Porter WR. 2012. Thermally accelerated degradation and storage temperature design space for liquid products. J Valid Technol 18(3):73– 92. 9. Li QC, Qiu F, Cohen K, Tougas T, Li J, McCaffrey J, Purdue T, Song JJ, Swanek F, Abelaira S. 2012. Best practices for drug substance stress and stability studies during early-stage development part I— Conducting drug substance solid stress to support phase Ia clinical trials. J Pharm Innov 7(3–4):214–224. 10. Porter WR. 2013. Degradation of pharmaceutical solids accelerated by changes in both relative humidity and temperature and combined storage temperature and storage relative humidity (T X h) design space for solid products. J Valid Technol 19(2). 11. Waterman KC, Swanson JT, Lippold BL. 2014. A scientific and statistical analysis of accelerated aging for pharmaceuticals, part 2: Assessment of error estimation methods. Manuscript submitted for publication. 12. Waterman KC, Swanson JT, Lippold BL. 2014. A scientific and statistical analysis of accelerated aging for pharmaceuticals, part 3: Shelf-life estimations when changes are within measurement error bars. Manuscript in preparation. 13. Porterfield RI, Capone JJ. 1984. Application of kinetic models and Arrhenius methods to product stability evaluation. Med Device Diag Ind 45–50. 14. King SP, Kung M, Fung H. 1984. Statistical prediction of drug stability based on non-linear parameter estimation. J Pharm Sci 73:657– 662. 15. Su XY, Po ALW, Yoshioka S. 1994. A Bayesian approach to Arrhenius prediction of shelf-life. Pharm Res 11(10):1462–1466. 16. Sundberg R. Statistical aspects on fitting the Arrhenius equation 1998. Chemometrics Int Lab Systems 41:249–252. 17. Gil-Alegre ME, Bernabeu JA, Camacho MA, Torres-Suarez AI. 2001. Statistical evaluation for stability studies under stress storage conditions. Il Farmaco 56:877–883. 18. Shimizu Y, Tamura T, Ono M, Kasai O, Nakajima T. 2002. Application of nonlinear fitting and selection of the most fitted equation by AIC in stability test of pharmaceutical ingredients. Drug Dev Ind Pharm 28(8):931–937. ´ cek L. 1997. Statistical properties off linearization 19. Kliˇca R, Kubaˇ of the Arrhenius equation via the logarithmic transformation. Chemometrics Intel Lab Sys 39:69–75. 20. Atkins PW. 1978. Physical chemistry. San Francisco: W.H. Freeman and Company, pp 868–869. 21. R Core Team. 2013. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. ISBN 3–900051–07–0. Accessed, at: http://www.R-project.org/ on 06 July 2013.
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
DOI 10.1002/jps.24075
A Scientific and Statistical Analysis of Accelerated Aging for Pharmaceuticals. Part 1: Accuracy of Fitting Methods KENNETH C. WATERMAN, JON T. SWANSON, BLAKE L. LIPPOLD FreeThink Technologies, Inc., East Lyme, Connecticut Received 18 February 2014; revised 31 May 2014; accepted 4 June 2014 Published online 8 July 2014 in Wiley Online Library (wileyonlinelibrary.com). DOI 10.1002/jps.24075 ABSTRACT: Three competing mathematical fitting models (a point-by-point estimation method, a linear fit method, and an isoconversion method) of chemical stability (related substance growth) when using high temperature data to predict room temperature shelf-life were employed in a detailed comparison. In each case, complex degradant formation behavior was analyzed by both exponential and linear forms of the Arrhenius equation. A hypothetical reaction was used where a drug (A) degrades to a primary degradant (B), which in turn degrades to a secondary degradation product (C). Calculated data with the fitting models were compared with the projected roomtemperature shelf-lives of B and C, using one to four time points (in addition to the origin) for each of three accelerated temperatures. Isoconversion methods were found to provide more accurate estimates of shelf-life at ambient conditions. Of the methods for estimating isoconversion, bracketing the specification limit at each condition produced the best estimates and was considerably more accurate than when extrapolation was required. Good estimates of isoconversion produced similar shelf-life estimates fitting either linear or nonlinear forms of the Arrhenius equation, whereas poor isoconversion estimates favored one method or the other depending on which condition C 2014 Wiley Periodicals, Inc. and the American Pharmacists Association J Pharm Sci 103:3000–3006, 2014 was most in error. Keywords: isoconversion; accelerated stability; shelf-life; Arrhenius fitting; kinetics; stability; chemical stability; mathematical model; in silico modeling
INTRODUCTION Determining the stability of drug substances and drug products is one of the critical, and at times, gating elements in the drug development process. In recent years, new protocols and mathematical models have been applied to the use of accelerated processes for rapidly determining the chemical stability of APIs and drug products both in solution and solid state.1–10 With the increased use of these models, it is appropriate to consider the implications of different methodologies in terms of both the statistics and science involved. In Part 1 of this detailed examination, the ability of different models to accurately estimate shelf-life is determined on a theoretical basis. Part 2 is focused on different statistical modeling methods for estimating the precision of projections.11 Finally, Part 3 examines the calculation options when there is little or no measured degradation.12 While the current paper examines what mathematical model best fits data without error bars, the second paper builds from that model to include statistical error propagation. The third paper evaluates the implications of different options within the model and error bar propagation of the first two papers for cases where there is little degradation occurring with respect to the error bars for the measurements. While previous examinations of accelerated stability modeling were based largely on loss of API potency with time,13–19 currently, increase in related substance (degradant) growth is more commonly the limiting factor on setting shelf-lives. This has implications in the modeling process for three reasons. First, the specification limit for the loss of potency (usually between 5% and 10% loss) is typically higher than the specification limit for Correspondence to: Kenneth C. Waterman (Telephone: +860-287-4253; E-mail: [email protected]) Journal of Pharmaceutical Sciences, Vol. 103, 3000–3006 (2014)
C 2014 Wiley Periodicals, Inc. and the American Pharmacists Association
3000
increase in the level of a degradant (typically 0.2%–1%). Second, the precision of the measurements for related substance growth is generally much tighter than that for potency loss. Third, in secondary degradation processes, potency loss can show simple kinetics even when individual degradant formation is complex. In classic chemical transformations, processes can be described by a reaction order related to the species involved in the rate-limiting step of the degradation process. An assumption is often made that the kinetics are pseudo-zero order because the extent of degradation is typically low (0.2%–1.0%) such that a straight line can be fitted to measured degradation versus time points. The rate under the measurement conditions will be the slope of that line; that is, the initial rate of degradant formation is often assumed to determine the rate constant. In solid-state and complex solution-state pharmaceutical systems, there are a number of reasons that this simple reaction order (and initial rate based) model does not describe the situation adequately, especially when focused on degradant formation rather than loss of API. First, degradants can undergo further degradation in secondary processes. This affects the primary degradant (rate slows as reaction proceeds) and secondary degradant (a lag followed by growth) formation kinetics. Second, diffusion processes can be rate-limiting steps, especially in solid-state chemical reactions. The kinetics of such reactions does not follow the simple reaction order model. Third, degradation products can catalyze further degradation in an autocatalytic process. This results in a lag followed by more rapid growth kinetics. A similar kinetic trend is seen in cases where an inhibitor (e.g., antioxidant) is consumed as the reaction proceeds. Finally, many solids are heterogeneous such that different nonequilibrating forms of the API react at different rates. Examples of this are drug molecules at interfaces versus in the bulk, and molecules at defect sites versus in crystalline domains. In our examination of calculation options, one of the evaluation factors for different
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
RESEARCH ARTICLE – Drug Discovery-Development Interface
methods will be the ability to accurately estimate behavior with complex kinetics. In experience with many solid drug products, we have seen that more than half do not follow a simple kinetic model. One of the processes for stability assessments evaluated here that has seen recent increased use in pharmaceutical applications is the isoconversion method.2–4 With isoconversion, rather than a focus on rate constants, the process involves an emphasis on the time hit the specification limit (the isoconversion time). For processes that can be characterized by a rate constant, the isoconversion time will be the specification limit divided by that rate constant. As a simple example for complex kinetics, we consider a process where a drug (A) goes to a primary degradation product (B), which undergoes further degradation to C, with the rate constants for the two processes corresponding to k1 and k2 : A→B→C
(1)
The amount of B and C as a function of time will correspond to Eqs. 2 and 3 (derived from Ref. 20): Bt = A0 Ct = A0 1 +
−k t k1 e 1 − e−k2 t k2 − k1
(2)
−k t 1 k2 e 1 − k1 e−k2 t k1 − k2
(3)
The behavior of B and C as a function of time is shown in Figure 1. As seen in the figure, B will rise with a decreasing rate, whereas C rises after a lag (as B builds up). Since the equations governing these species are known, one can compare the results from an exact solution of the problem to those obtained using the different methods for using accelerated data to estimate ambient stability. The kinetic forms reflected in Eqs. 2 and 3
3001
will be representative of many of the types of observed non-zero order behaviors in drug products. The dependence of solid-state chemical reactions on temperature and RH has been shown to follow a modified Arrhenius equation which can be expressed in linear (Eq. (4) and exponential (Eq. (5) forms. ln k = ln A − Ea /(RT) + B(RH)
(4)
k = A e(−Ea /RT)+B(RH)
(5)
where k is a rate for degradant formation, T is the temperature (K), and RH is the equilibrium relative humidity to which the sample is exposed (as opposed to the external RH for a packaged product).1,2 Varying T and RH to determine the rate at each condition, then determining the best fit of the variables A, Ea , and B allows a mathematical model to be built which can be used to calculate the rate (and corresponding shelf-life) at extrapolated or interpolated T and RH conditions. In real applications of accelerated aging, the RH dependence complicates predictions of ambient behavior yet cannot be ignored. Because the form of the equations remains the same with or without the humidity term, the calculations herein assume a constant RH throughout with the implication that the preferred methodologies for accurate predictions of shelf-lives will be the same when the RH term is added. The added complexity of the modified Arrhenius equations, however, will exacerbate any inaccuracies intrinsic in a method. For purposes of illustration, we are using a system where the true kinetic form of the degradant formation is known. In practice, this is often not the case. One method for determining shelf-life when the kinetic form is not known involves calculation of the isoconversion time; that is, the time when the degradant reaches its specification limit at each condition. At the specification limit, the isoconversion time will be proportional to the (isoconversion) rate. This allows the isoconversion time to be used to calculate the Arrhenius parameters at the specification limit, without knowing the form of the rate equation. The shelf-life at ambient conditions is the isoconversion time at those corresponding conditions. One can illustrate the relationship mathematically for several simple cases. In measuring degradant formation (d) for a zero-order or pseudo-zero order reaction with a rate constant k, degradant formation will be proportional to time (t) as shown in Eq. (6): d = −kt
(6)
The time at which the degradant reaches the specification limit is defined as the isoconversion time (iso). Since d at each isoconversion time is the same at all temperatures (it is the specification limit), the following will be true: kT1 iso1 = kT2 iso2 Figure 1. Model drug behavior used in the present method evaluation showing the 60◦ C curves. Drug substance A degrades to a primary degradant B (blue curve), which in turn degrades to a secondary degradant C (orange curve). The graph shows B and C as a function of time and typical specification limits of 0.2% and 0.5% for degradant formation. DOI 10.1002/jps.24075
(7)
Rearranging this equation gives the result that the ratio of isoconversion times is inversely proportional to the rate constants: kT iso2 = 1 iso1 kT2
(8)
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
3002
RESEARCH ARTICLE – Drug Discovery-Development Interface
One can substitute into the rate equation (Eq. (6) the specification limit, the rate constants (as determined by the Arrhenius equation) and isoconversion times at two temperatures to determine the activation energy (Eqs. 9–12). spec = A e−Ea /RT1 iso1
(9)
spec = A e−Ea /RT2 iso2
(10)
E iso2 − a =e R iso1
1 T1
R ln Ea = −
1 T1
− T1
iso2 iso1
−
2
1 T2
(11)
(12)
Once the activation energy is determined, one can calculate the ratio of rates between a known temperature and any other temperature. Determining the isoconversion times at two temperatures is therefore the minimum requirement to determine the shelf-life at ambient conditions. The true rate does not need to be known. This can be extended to more complex kinetics. For example, first order kinetics follows Eq. (13): d = [A]0 1 − e−kt
(13)
d − [A]0 = −e−kt [A]0
(14)
or
Taking the natural logarithm gives Eq. (15). − ln
d − [A]0 [A]0
= −kt
(15)
This is the same as Eq.
(6) with d replaced by the expresd−[A]0 0 , and − ln will be the same at differsion − ln d−[A] [A] [A] 0
0
ent temperatures at the specification limit. Therefore, Eq. (12) again allows one to relate the isoconversion time at an accelerated temperature to the isoconversion time at an ambient condition (shelf-life). In fact, no matter what the form of the degradant formation kinetics, using isoconversion will maintain this relationship, assuming the reaction mechanism to form the degradant does not change with temperature. This means that the shelf-life can be extrapolated from accelerated conditions without knowing the true kinetic form for degradant formation as long as the isoconversion times are known at two or more temperatures. For the present evaluation, rather than illustrate a more complex example using equations, several methods of shelf-life prediction will be compared with the value one would obtain if the true kinetic form was known. Three methods for estimating degradant formation rates are compared with respect to their ability to provide accurate estimations of shelf-life from accelerated data.
First, the amount of degradant at each time point is considered separately with respect to the initial amount of degradant. Effectively, each measured value at time ti and temperature Ti will result in a rate constant ki based on the change in the amount of degradant divided by time. These values can be used to determine a set of independent rates based on all time points and conditions. Second, each temperature condition is used to determine a best fit, zero-order rate constant for that condition. This is effectively a rate based on the initial slope for degradant formation. Multiple time points contribute to the rate estimation in each condition. Third, an isoconversion method is employed. In practice, the isoconversion time point will not be known a priori, and multiple measurements of degradant formation as a function of time will be made. These measurements can be used to estimate the isoconversion time at each condition. How best to use these data for the estimation of isoconversion will be examined using three methods. (1) A line will be drawn through the data to find the intercept with the specification limit. (2) The two closest points to the specification limit will be used to either interpolate or extrapolate to the specification limit. (3) The single point closest to the specification limit will be used to define a rate using k = d/t. Each of the methods results in a different rate or set of rates at a given condition, which will have implications when these rates are used with different temperatures to predict shelflife, particularly when the extent of reaction is different at the different temperatures. In those cases where ideal linear (zero order) behavior, without noise, is observed (or when the rate can be characterized effectively based on the initial slope), all the methods will result in the same rate constant. In our examination of drug product stability, over half the cases studied (n > 200) do not show such ideal behavior. In addition, since most pharmaceutical stability studies used in product development involve very few time points, it can be difficult to validate linearity in a given system. We will examine cases where there is apparent linearity, but where the degradant level remains well below the specification limit to compare how accurately the various treatments predict ambient behavior in such cases. Two methods are examined for determining the temperature dependence of degradant formation and thereby estimating the shelf-life at ambient conditions. Each assumes the temperature dependence can be modeled by the Arrhenius equation. (1) The equation can be transformed to a linear form (Eq. (4) and parameters determined by a linear least squares method. (2) The equation (Eq. (5) can be used directly and parameters determined using a nonlinear least-squares method. In either case, one is minimizing an objective function, typically the sum of the squared error between the actual and predicted values. If the values fit the theoretical model perfectly, both methods will give the same results. With real data, however, this situation will not be true, and the two methods will produce different results because of the difference in the objective function. When fitting Arrhenius data, the relative leverage of the higher temperature points will be greater for the nonlinear fit, since the predicted value depends exponentially on the temperature. Loss of potency can be particularly difficult in some cases for Arrhenius fitting since multiple pathways for loss of API can each have different activation parameters; however, this issue is not evaluated in this study. Similarly, differences in the
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
DOI 10.1002/jps.24075
3003
RESEARCH ARTICLE – Drug Discovery-Development Interface
Table 1.
Primary Degradant-Based Shelf-Life with Different Modelsa Shelf-Life (yrs) 25◦ C
Method Specification Limit 0.2%
Exactb Single constant time point (14 days)c Individual time pointsd Line through all pointse Single point at isoconversionf Linear, all-point isoconversiong Bracketed isoconversionh Closest point isoconversioni
Specification Limit 0.5%
Linear
Nonlinear
Linear
NonLinear
1.43 0.56 0.62 0.29 1.43 12.35 1.36 1.23
1.43 0.40 0.54 0.19 1.43 467.94 1.26 1.25
4.45 1.40 1.56 0.71 4.45 1.40 3.19 2.76
4.45 1.00 1.36 0.49 4.44 1.16 4.39 2.63
a Shelf-life at 25◦ C estimated from accelerated conditions of 50◦ C, 60◦ C, and 70◦ C for a hypothetical drug product with the shelf-life being limited by formation of a primary degradant that can degrade to a secondary degradation product (k1 and k2 at 50◦ C equal 0.000113 and 0.011250%/d). Time points at each temperature, in addition to the origin, were 3, 7, 14, and 28 days b Refers to the true value based on the rate constants. c Based on linear rate constants using only one time point at each temperature. d Based on four linear rate constants at each temperature. e Based on using linear rate constants through four time points at each temperature. f Using the exact time points corresponding to the intersection of the degradant formation curves with the specification limit. g Estimating isoconversion times by linear fitting of all four time points at each condition to determine the intersection with the specification limit. h Estimating isoconversion times using the two time points closest to the specification limit at each temperature (or extrapolating from the last two time points when the specification limit is greater than the values at that temperature). i Estimating isoconversion times using the point (one of four) closest to the specification limit at each temperature.
activation parameters for the two rate constants in our A to B to C scenario are not considered here. In the present investigation, the focus is on the impact of these different calculation methods on the accuracy of predictions using perfect (i.e., no error) data, but with complex (but common) kinetics. It can be asserted that the ability of a method to provide accurate predictions for these types of kinetics will extend to other complex (or simple) kinetic schemes. We assume that the true kinetic form is unknown in the investigation, which is generally the case in practice, since it is rare that stability studies are carried out with enough time points and duration to truly determine the kinetic form. Even when kinetics appears linear over a range, curvature at greater conversions (nearer to the specification limit) may significantly affect the estimations.
METHODS Many of the calculations were carried out using the R statistical package.21 The R statistical package contains a number of subroutines used herein, including the “lm function,” which is a linear fitting module, and “nls function,” which is a nonlinear, least squares fitting routine. Both these routines are available as part of this open-source package. Values of the percent degradant as a function of time were calculated using Eqs. 2 and 3. The rate constants used met the following three criteria: the isoconversion time with an 0.2% specification limit and a temperature of 50◦ C is within the 28 day window of the hypothetical accelerated stability study; the isoconversion time with an 0.5% specification limit and a temperature of 50◦ C exceeds that window (approximately 60 days); the shelf-life at 25◦ C with a 0.5% specification limit is approximately 4 years (activation energy of 25 kcal/mol). Rates at 60◦ C and 70◦ C were determined assuming the 25 kcal/mol activation energy. For the calculations of B in the A to B to C scheme, k1 and k2 at 50◦ C equal 0.000113 and 0.011250%/d, respectively; for calculations of C, k1 , and k2 at 50◦ C equal 0.000112 and 0.090000%/d, DOI 10.1002/jps.24075
respectively. The points used for the calculations at each temperature were 0, 3, 7, 14, and 28 days. Three methods were used to calculate the rate constants used for determining the Arrhenius parameters. The first method calculates the rate of each individual time point as k = d/t. This method was used for calculations with a single constant time point (d was the value determined at a time of 14 days), a single measured point closest to isoconversion, and multiple individual time points. The second method fit the points to a line using the lm function of the R statistical package.21 The third method calculated k as k = specification limit/isoconversion time. This method was used for the single point at isoconversion, the linear fit to isoconversion, and the bracketed isoconversion. In the latter two cases, the isoconversion time was first calculated by determining the intercept of the best fit line with the line d = specification limit. The lm function was used to fit the points (all or the bracketing points). The Arrhenius parameters were determined by fitting Eq. (16) (linear fit) using the lm function or Eq. (17) (nonlinear fit) using the nls function of the R statistical software package. ln k = ln A −
Ea RT
k = Ae−Ea /RT
(16)
(17)
The Arrhenius parameters were used to calculate k at 25◦ C, which was then used to calculate the shelf-life as specification limit/k.
RESULTS AND DISCUSSION Table 1 summarizes the results of the calculations using different methods for determining the extrapolated shelf-life of the B degradant in an A to B to C scheme. Results reflect the use of Eq. (2) to calculate degradant formation. The shelf-life was
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
3004
RESEARCH ARTICLE – Drug Discovery-Development Interface
Table 2.
Secondary Degradant-Based Shelf-Life with Different Modelsa Shelf-life (yrs) 25◦ C
Method Specification Limit 0.2%
Exactb Single constant time point (14 days)c Individual time pointsd Line through all pointse Single point at isoconversionf Linear, all-point isoconversiong Bracketed isoconversionh Closest point isoconversioni
Specification Limit 0.5%
Linear
Nonlinear
Linear
Nonlinear
2.02 8.29 16.64 3.29 2.02 2.75 2.06 2.14
2.02 4.24 6.84 2.15 2.02 2.05 2.05 3.12
4.01 20.74 41.61 8.21 4.01 7.56 4.78 7.41
4.01 10.60 17.11 5.36 4.01 5.23 4.25 6.56
a Shelf-life at 25◦ C estimated from accelerated conditions of 50◦ C, 60◦ C, and 70◦ C for a hypothetical drug product with the shelf-life being limited by the formation of a secondary degradant (k1 and k2 at 50◦ C equal 0.000112 and 0.090000%/d). Time points at each temperature, in addition to the origin, were 3, 7, 14, and 28 days b Refers to the true value based on the rate constants. c Based on linear rate constants using only one time point at each temperature. d Based on four linear rate constants at each temperature. e Based on using linear rate constants through four time points at each temperature. f Using the exact time points corresponding to the intersection of the degradant formation curves with the specification limit. g Estimating isoconversion times by linear fitting of all four time points at each condition to determine the intersection with the specification limit. h Estimating isoconversion times using the two time points closest to the specification limit at each temperature (or extrapolating from the last two time points when the specification limit is greater than the values at that temperature). i Estimating isoconversion times using the point (one of four) closest to the specification limit at each temperature.
then calculated for two different specification limits, with different modeling methods using both the linear and nonlinear versions of the Arrhenius equation. Similarly, Table 2 provides the results for the C degradant. Formation of a primary degradant that undergoes further degradation to a secondary degradant is a common example of complex degradation kinetics. Traditionally, stability projections using accelerated data have assumed linearity (zero order or pseudo zero order kinetics, or use of initial slope). As seen in the results for the primary degradant (Table 1), whether using a single time point or multiple points, this assumption leads to significantly inaccurate estimations of room temperature behavior. With formation of a primary degradant, the linearity assumption underestimates the ambient shelf-life. This occurs because the lower temperature conditions do not show the downward curvature to as great an extent as with the higher temperatures. This leads to a lower apparent temperature slope (activation energy) and a rate higher than the true rate at ambient conditions. With formation of a secondary degradant, there is a lag before growth. The results with an assumption of linearity (Table 2) are again quite inaccurate. In this case, the lag phase is not yet exhausted at lower temperatures. Since higher temperatures show greater reactivity, the linear approximation results in an apparent higher activation energy and an overestimation of the shelf-life at ambient conditions. While the primary degradant would give an estimation that is overly conservative, with a secondary degradant, the linear approximation could result in product failures on storage. As long as the degradant formation (whether primary or secondary) is assumed to be linear (zero order or pseudo zero order), the estimations to ambient conditions from high temperatures are not in general improved by using a nonlinear least squares fitting to the Arrhenius equation. The rates themselves at the accelerated conditions do not accurately represent the behavior to the specification limit. The form of the Arrhenius equation that is fit cannot fix this problem.
The situation is significantly improved when using isoconversion times. With the exact isoconversion times used for each condition, the projections to ambient stability are exactly correct for both primary and secondary degradants (Tables 1 and 2) as anticipated in the introduction. The practical challenge for using accelerated aging in such systems is how to determine the isoconversion times when they are not known a priori. Three different approaches at two specification limits were evaluated in this study. There are two reasons for examining multiple specification limits. Whereas with the zero-order assumption one has the same rate at all times, this is not true for complex kinetics. Because the specification limit is the time at which the degradant reaches the specification limit, the effective rate will depend on the specification limit and therefore the Arrhenius parameters (A, Ea ) will also depend on the specification limit. The two cases also illustrate another important difference. At the lower specification limit, all isoconversion times can be interpolated. At the higher specification limit, extrapolation is needed for the 50◦ C point. This allows comparison of the accuracy of interpolated and extrapolated data. The most obvious method for estimating isoconversion is to use a linear fit of the time versus percent degradant data to determine the time to reach the specification limit. The isoconversion rate would then be the difference between the specification limit and the initial degradant level divided by the isoconversion time. This is effectively the same as the linear kinetics assumption (zero or pseudo zero order), except that here the rate propagated to the Arrhenius fitting is the rate to the specification limit. This method creates a significant error when the fitted line is a poor estimation of the true isoconversion. This will happen when most of the points used to fit the line are well beyond isoconversion or extrapolation is required to reach the isoconversion time. Whether the isoconversion time is overestimated or underestimated depends on the curvature of the degradant time plot and whether the isoconversion time is being interpolated or extrapolated. These factors also influence
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
DOI 10.1002/jps.24075
RESEARCH ARTICLE – Drug Discovery-Development Interface
Figure 2. When including points well past isoconversion (illustrated here with degradant B at 70◦ C), a linear fit to the data may give a nonzero intercept. In this example, this results in underestimating the isoconversion time (estimated isoconversion time of 0.56 days shown in red versus the true isoconversion time of 2.04 days, shown in green).
whether it is the higher or lower temperature points which have the least accurate estimate of the isoconversion time. A second strategy to estimate isoconversion evaluated here is to use only the measured points that bracket the specification limit (or extend the closest two for extrapolations), rather than use all the points measured. This better accounts for the fact that the rate varies with time. As can be seen in Tables 1 and 2, the estimations with this process are significantly better for both the primary and secondary degradants. While not perfect, the “bracketed” approach provides a reasonable estimation of the shelf-lives in all cases. Comparing the predictions from the different methods to estimate isoconversion in the examples illustrates several key points. First, using all four time points works particularly poorly when estimating the shelf-life of the primary degradant at a 0.2% specification limit. This is the result of a single poor estimation—the isoconversion time at 70◦ C (Fig. 2). When points well past the isoconversion time are used, as in this ex-
3005
ample, the isoconversion time is underestimated and a serious overestimation of the shelf-life results. This problem is easily remedied by not including points significantly past the isoconversion time. A more difficult problem occurs when the study is not carried to isoconversion, as happens at 50◦ C with the 0.5% specification limit (Fig. 3). In this example, even though the points measured show nearly perfect linearity (R2 = 0.998), the nonlinearity of the true kinetic form is observed by the time isoconversion is reached. Linearity observed over the range of the experiment is insufficient to assume zero-order kinetics. Degradant formation must remain linear out to the isoconversion time. The result is a significant underestimation of the isoconversion time at 50◦ C, which results in a corresponding underestimation of the shelf-life. In general, when kinetics is complex, extrapolated values can be much less accurate than interpolated values, which are reflected in the estimate of the shelf-life. Bracketing the isoconversion time works well in these examples because the data are exact. In real studies, one must balance minimizing the error that results with complex kinetics when a linear estimation is used, with the error inherent in measurement of the data points. In the former case, focusing on the isoconversion time will reduce the error. In the latter case, more measurements will provide a more accurate estimate of the true degradant levels as a function of time. The commercial software ASAPprime incorporates an algorithm to look at how many points are appropriate to use for the estimation in the event the overall fit cannot be approximated as linear. The final estimation method for isoconversion uses just the point closest to the specification limit to define a rate. As seen in Tables 1 and 2, the estimations from this method are not as good as with the “bracketed” process; however, they are better than most other methods. The fact that a single point near isoconversion provides a better estimation than most of the other methods illustrates the importance of accurately estimating the isoconversion time for complex kinetics. This is not a method we would recommend, as per the preceding discussion, it benefits in this simple example from the fact the data points exactly fit the kinetic scheme. R
Figure 3. When measurements are not carried out to isoconversion (illustrated here with degradant B at 50◦ C), curvature may not be apparent (R2 for the fitted line is 0.998). In this example, the isoconversion determined by a linear fit (estimated isoconversion time of 51.6 days shown in red versus the true isoconversion time of 62.0 days, shown in green), emphasizing the desirability of carrying out stability studies to the specification limit at each condition. DOI 10.1002/jps.24075
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
3006
RESEARCH ARTICLE – Drug Discovery-Development Interface
When nonlinear fitting to the Arrhenius equation is employed with the isoconversion process, the results in most cases are not significantly different from those with the linear fitting. When the estimation method is exactly correct, the linear and nonlinear fits produce the same result. When the estimation method is not exact, since the nonlinear fit leverages the higher temperature points more heavily than the linear fit, it will tend to give better estimations when the high temperature estimations of isoconversion are closer to the true isoconversion values. The opposite will be the case when the low temperature values are better estimations. In the examples used, the higher temperature conditions required only interpolation to determine the isoconversion times. The lower temperature conditions in some cases required extrapolation. The result is that in some cases the nonlinear fit appears to work better. The more significant difference between these methods relates to the estimated uncertainties of prediction produced by these methods, which is outside the scope of the present paper (discussed in the subsequent papers11,12 ). It is possible to include a weighting factor onto either linear or nonlinear least-squares fitting methods. In particular, since nonlinear (exponential) fitting biases toward the high temperature points, it should be possible to counter this by weighting the lower temperature points more heavily. Ultimately, this opposes the argument that using the “true” form of the Arrhenius equation intrinsically should be more accurate.14 While it is likely that with appropriate weighting nonlinear methods could approximate the linear methods; however, since the latter are computationally less complex, it makes more sense to simply use linear fitting. Moreover, any weighting will also play a part in the statistics for the uncertainty of the prediction, which is discussed in the subsequent paper in this series.11
CONCLUSIONS Traditional (i.e., linear rate constants estimated by initial slope) and isoconversion methods were compared for the accuracy of predicted shelf-lives from accelerated data when the kinetics of primary and secondary degradant formation was complex. The results show clearly that traditional stability assumptions provide poor estimates. This is true whether a single or multiple time points are used, and whether a linear or nonlinear Arrhenius fitting is used. With isoconversion calculations, shelf-life predictions are directly related to how good an estimate of isoconversion is achieved at each condition. Of the methods examined, bracketing of the specification limit (extending the last two time points for extrapolations) provides the best predictions. The best results are obtained when the isoconversion estimates are obtained by interpolation. Even relatively short extrapolations can introduce considerable error. In this study, linear and nonlinear fitting (parameter estimation) of the Arrhenius equation provide similar accuracy for ambient shelf-life predictions. A key point is that if the isoconversion times are accurate, the two approaches produce the same result. It is not possible to decide whether linear or nonlinear fitting provides better accuracy since the results depend on which condition is in error. A proper discussion of these two approaches requires also considering the estimate of uncertainty they produce (as discussed in the second paper in this series11 ). Conflict of interest: The authors have a commercial stake in the software ASAPprime , which broadly uses accelerated aging data with some of the methods described in the manuscript. R
REFERENCES 1. Waterman KC, Adami RC. 2005. Accelerated aging: Prediction of chemical stability of pharmaceuticals. Inter J Pharm 293:101–125. 2. Waterman KC, Carella AJ, Gumkowski MJ, Lukulay P, MacDonald BC, Roy MC, Shamblin SL. 2007. Improved protocol and data analysis for accelerated shelf-life estimation of solid dosage forms. Pharm Res 24(4):780–790. 3. Waterman KC, Colgan ST. 2008. A science-based approach to setting expiry dating for solid drug products. Regul Rapporteur 5(7):9–14. 4. Waterman KC. 2009. Understanding and predicting pharmaceutical product shelf-life. In: Huynh-Ba K; Ed. Handbook of stability testing in pharmaceutical development. New York: Springer, pp 115– 135. 5. Waterman KC, MacDonald BC. 2010. Package selection for moisture protection for solid, oral drug products. J Pharm Sci 99:4437–4452. 6. Waterman KC. 2011. The application of the accelerated stability assessment program (ASAP) to quality by design (QbD) for drug product stability. AAPS PharmSciTech 12:932–937. 7. Colgan ST, Watson TJ, Whipple RD, Nosal R, Beaman JV, De Antonis DM. 2012. The application of science- and risk-based concepts to drug substance stability strategies. J Pharm Innov 7(3–4):205–213. 8. Porter WR. 2012. Thermally accelerated degradation and storage temperature design space for liquid products. J Valid Technol 18(3):73– 92. 9. Li QC, Qiu F, Cohen K, Tougas T, Li J, McCaffrey J, Purdue T, Song JJ, Swanek F, Abelaira S. 2012. Best practices for drug substance stress and stability studies during early-stage development part I— Conducting drug substance solid stress to support phase Ia clinical trials. J Pharm Innov 7(3–4):214–224. 10. Porter WR. 2013. Degradation of pharmaceutical solids accelerated by changes in both relative humidity and temperature and combined storage temperature and storage relative humidity (T X h) design space for solid products. J Valid Technol 19(2). 11. Waterman KC, Swanson JT, Lippold BL. 2014. A scientific and statistical analysis of accelerated aging for pharmaceuticals, part 2: Assessment of error estimation methods. Manuscript submitted for publication. 12. Waterman KC, Swanson JT, Lippold BL. 2014. A scientific and statistical analysis of accelerated aging for pharmaceuticals, part 3: Shelf-life estimations when changes are within measurement error bars. Manuscript in preparation. 13. Porterfield RI, Capone JJ. 1984. Application of kinetic models and Arrhenius methods to product stability evaluation. Med Device Diag Ind 45–50. 14. King SP, Kung M, Fung H. 1984. Statistical prediction of drug stability based on non-linear parameter estimation. J Pharm Sci 73:657– 662. 15. Su XY, Po ALW, Yoshioka S. 1994. A Bayesian approach to Arrhenius prediction of shelf-life. Pharm Res 11(10):1462–1466. 16. Sundberg R. Statistical aspects on fitting the Arrhenius equation 1998. Chemometrics Int Lab Systems 41:249–252. 17. Gil-Alegre ME, Bernabeu JA, Camacho MA, Torres-Suarez AI. 2001. Statistical evaluation for stability studies under stress storage conditions. Il Farmaco 56:877–883. 18. Shimizu Y, Tamura T, Ono M, Kasai O, Nakajima T. 2002. Application of nonlinear fitting and selection of the most fitted equation by AIC in stability test of pharmaceutical ingredients. Drug Dev Ind Pharm 28(8):931–937. ´ cek L. 1997. Statistical properties off linearization 19. Kliˇca R, Kubaˇ of the Arrhenius equation via the logarithmic transformation. Chemometrics Intel Lab Sys 39:69–75. 20. Atkins PW. 1978. Physical chemistry. San Francisco: W.H. Freeman and Company, pp 868–869. 21. R Core Team. 2013. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. ISBN 3–900051–07–0. Accessed, at: http://www.R-project.org/ on 06 July 2013.
Waterman, Swanson, and Lippold, JOURNAL OF PHARMACEUTICAL SCIENCES 103:3000–3006, 2014
DOI 10.1002/jps.24075