- Email: [email protected]
A mathematical model describing the thermal virus inactivation
Vaccine 19 (2001) 3575– 3582 www.elsevier.com/locate/vaccine
A mathematical model describing the thermal virus inactivation Catherine Noe¨l *, Sandrine Charles, Alain Franc¸on, Jean Pierre Flandrois De6elopment Department, A6entis Pasteur, 69280 Marcy l’Etoile, France Received 1 May 2000; received in revised form 4 December 2000; accepted 12 December 2000
Abstract A new mathematical model is proposed to describe the inactivation of viruses at different temperatures. This model takes into account the exponential decrease of the viral titer with time, the inactivation rate being an exponential function of the temperature. A one-step non-linear regression was used to fit oral poliovirus vaccine (OPV) experimental data. In one of the applications of the model, we illustrate the use of our model to compare the accelerated degradation test of OPV new formulations to standard OPV. Such a model is both simple and convenient to use. It should be a useful tool in optimizing formulations for live viral vaccines. © 2001 Elsevier Science Ltd. All rights reserved. Keywords: Kinetic; Mathematical model; Oral polymyelitis vaccine; Temperature; Thermal inactivation; Viral vaccines
1. Introduction Stability of biological products is a major requirement in the pharmaceutical industry [1], ensuring quality, safety of use and, for vaccines success of the immunization. Temperature is the major environmental factor that affects vaccine stability [2]. It is studied because of its practical importance as a limiting factor of the stability of live viral vaccine against polio (OPV), measles, mumps, rubella, varicella. The common method used to study the thermostability of live virus vaccines is to analyze the loss of titer (i.e. the decrease in virus number) with time, during storage at elevated temperatures for which inactivation will occur more rapidly. This is the so-called accelerated inactivation test or accelerated degradation test. The test usually consists of two steps (i) description of the exponential decrease of the titer with time [3–9], by computing the inactivation rate versus time after linearization; and (ii) the use of the Arrhenius equation which correlates the natural logarithm of the inactivation rate with the inverse of temperature [7,10 – 12].
* Corresponding author. Tel.: + 33-4-37373456; fax: + 33-437373029. E-mail address: [email protected] (C. Noe¨l).
This two-step approach is up to now not satisfactory; the error associated with the determination of viral titer are not taken into account, and therefore, the confidence interval parameters may be too wide to apply the model in a reliable predictive manner [13,14]. A more sophisticated way is to establish a mathematical relationship to directly describe the decrease of the viral titer as a function of both time and temperature during storage. We previously tried to build such a model based on the Arrhenius equation under its exponential form but it was impossible to fit the model to our data sets. The Arrhenius equation is based on a determinist mechanistic approach. On the other hand, we propose in this work an empirical approach with a formulation of the model based on simple mathematical functions, chosen for their simplicity and their suitability to fit our data. Our aim was to accurately describe the biological phenomenon with a robust and simple model, which can be used, for example, to compare new formulations of live viral vaccines. We, thus, developed this model on oral polioviruse vaccine and we applied it to compare the stability of different formulations. This paper focuses first on the non-linear fit of the mathematical function to experimental data (either from bioassays or literature), and second on the use of the model in comparison with virus inactivation in various formulations.
0264-410X/01/$ - see front matter © 2001 Elsevier Science Ltd. All rights reserved. PII: S0264-410X(01)00011-1
C. Noe¨l et al. / Vaccine 19 (2001) 3575–3582
3576
2. Materials and methods
2.1. Data We studied the inactivation of the three strains of Polio Sabin virus included in the oral polio vaccine (OPV), all of which are stored in liquid form. The data were either taken from the literature or obtained in Aventis Pasteur laboratories. Data set 1 was used to fit the model. It was extracted from Husson et al. [15]. The authors describe a trivalent
Table 1 Sampling schedule for Data set 2 and 3 to study the inactivation of polioviruses Data set
T (°C)
Sampling time
2
5 25 37 45 5 25 37
1, 3, 6, 9, 12, 18, 24, 30, 36 months 2, 4, 6, 8, 10, 12 weeks 2, 5, 7, 14, 21 days 2, 5, and 7 days 1, 3, 5, 6, 7 months 2, 4, 6, 8, 10, 12, 16 weeks 5, 7, 9, 10, 11, 12, 14, 16, 18, 19, 20, 21 days 2, 4, 6, 8, 10 days 1, 2, 3, 4, 5, 6 days 0.62, 0.83, 1.06, 1.25, 1.67 days
3
40 45 50
working reference preparation of OPV containing the three Sabin Strains. The global titer was measured after storage at 20°C for 1, 2, and 3 weeks; at 37°C for 1, 3 and 5 days; and at 45°C for 1, 2 and 3 days. The titer given in the table of that paper is the mean of two measures. An estimation of the titer at the beginning of the study time (t= 0) is given by the mean of 68 measures of the reference stored at − 70°C. Data sets 2 and 3 were used to fit the model. They were obtained in Aventis Pasteur laboratories. They provided the attenuated poliovirus strains used in this investigation: type 1: LSC 2ab/KP2 So + 2; type2: P712 ch 2ab/KP3 So +2; type3: Leon 12a1b/KP4 SO +4. In Data set 2, the titer was studied for OPV exposed to 5, 25, 37 or 45°C for 1–36 months (Table 1). Data set 3 concerns OPV exposed to six temperatures ranging for 5–50°C for 1–7 months (Table 1). Viral titer was defined in CCID50 per dose of vaccine. The titer of each strain was measured after neutralization of the two other strains using the appropriate antisera. Serial dilutions of virus suspension were performed using the culture medium, each dilution was inoculated in ten wells of a microplate and HepII cells were added in each well. The microplates were incubated at 32°C with CO2 and the cytopathogenic effect observed to determine the virus titer. Calculations were done using the Spearman–Karber method.
2.2. Comparison data set Data sets 4 used for comparison of various formulations were extracted from Leal et al. [12]. Different formulations of trivalent OPV were tested, each constituting a data set. The titer was analyzed after 4, 7, 14, 21, 28, 35 and 42 days at 20°C. 4, 7, 14, 21, 28 and 35 days at 26°C. 1, 2, 4, 7, 14 and 21 days at 37°C. The three virus strains were studied. Titers given in tables of the article represent the mean of four measures.
3. Modelling
3.1. Model building The inactivation of biological product exposed to heat is usually described by the following model [3–9]: Fig. 1. Filled contour plot of the global model Eq. (6) describing the effect of (h,i) on the titer loss (lnN0 − lnN), for a temperature of 20°C and a sampling time of 19 days: the titer losses are increasing from the black area to the white area. For example if h= 0.0005 and i = 0.05 the estimated titer loss is 0.02 log DICC50/ml and if h = 0.002 and i = 0.20 the estimated titer loss is 2.07 DICC50/ml when the vaccine is stored 19 days at 20°C. Moreover, all the points placed on the same line had an equal predicted loss.
Nt = N0 exp(− kt) U log Nt = log N0 − kt U ln Nt =ln N0 − kt
(1)
where Nt refers to the titer at time t, N0 to the initial titer, and k to the inactivation rate for a given temperature T. k has a biological interpretation as the inactivation rate (i.e. the loss of titer per unit of time).
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
3577
Fig. 2. Fit of the new global model Eq. (3) with Data set 2. Squares represent observed values, the continue-line the adjusted global model and the dash lines the 95% confidence limits calculated with Eq. (5).
Fig. 3. (a) Standardized residuals versus the calculated values of Data set 2; such a plot corroborates the good fit of the model; no obvious heteroscedasticity can be observed and the residuals seem to be randomly distributed. (b) Normal quantiles versus standardized residual quantiles of Data set 2; this plot shows the normality of the residual distribution and the absence of outliers.
To build the model, we studied the relation between k and the temperature T. We then propose a new relationship between k and T expressed as follows: k = h exp(iT)
(2)
where h correspond to the inactivation rate for temperature equal to 0, and i have no biological signification.
By substitution of Eq. (2) into Eq. (1), we finally get a complete model. ln N(t,T) = f(t, T, U) ln N(t,T) = ln N0 − h exp(iT)t
(3)
or ln N0 − ln N(t, T) = h exp(iT)t
(4)
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
3578
Table 2 Asymptotic results obtained by fitting model IV to Data sets 1, 2 and 3: parameter estimations, their 95% standard error, the t ratio and the P-value, and the correlation matrixa Fit data set
Virus titration
Parameter
Estimate
Standard Error
1 (n= 12)
Global
ln N0 h i
12.09 0.00800 0.127
ln N0 h i
2 (n = 81)
2 (n= 81)
2 (n=81)
3 (n= 146)
3 (n=134)
3 (n=149)
a
Strain: 1
Strain:2
Strain: 3
Strain: 1
Strain: 2
Strain: 3
t ratio
P value
Correlation matrix
Iteration number
0.26 0.00333 0.009
46.4 −2.4 13.9
B0.01 0.04 B0.01
Á 1 Ã Ã 0.58 Ã Ä−0.47
1 −0.98
 à à à 1Å
18
14.45 0.00080 0.149
0.10 0.00011 0.003
152.1 −7.4 47.3
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.61 Ã Ä−0.41
1 −0.93
 à à à 1Å
6
ln N0 h i
12.41 0.00065 0.150
0.10 0.00011 0.004
123.5 −5.9 37.7
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.61 Ã Ä−0.42
1 −0.94
 à à à 1Å
5
ln N0 h i
14.19 0.00085 0.141
0.10 0.00013 0.004
138.8 −6.8 39.2
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.60 Ã Ä−0.31
1 −0.89
 à à à 1Å
7
ln N0 h i
13.75 0.00083 0.163
0.11 0.00008 0.002
127.1 −10.3 69.1
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.31 Ã Ä−0.08
1 −0.96
 à à à 1Å
8
ln N0 h i
10.70 0.00118 0.154
0.14 0.00014 0.003
78.2 −8.7 56.2
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.35 Ã Ä−0.07
1 −0.94
 à à à 1Å
8
ln N0 h i
13.40 0.00107 0.153
0.15 0.00015 0.004
91.4 −7.0 43.0
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.27 Ã Ä−0.02
1 −0.96
 à à à 1Å
7
The iteration number needed to estimate various parameters is also indicated.
Table 3 Asymptotic results obtained by fitting model IV to data sets 4: parameter estimations and their 95% standard error h Calculated
lnN0
i Calculated
Calculated
Observed
TYPE 1 RIT (reference = 1) SU-ARG (2) MgCl2-ARG (3) MgCl2-199 (4)
13.71 13.83 13.80 13.96
[13.83;13.95] [13.70; 13.95] [13.62; 13.97] [13.79; 14.12]
13.91 13.84 13.86 14.02
0,00079 0.00180 0.00133 0.00196
[0.00044; [0.00126; [0.00073; [0.00122;
0.00113] 0.00234] 0.00192] 0.00270]
0.166 0.149 0.159 0.146
[0.155; [0.140; [0.147; [0.135;
0.178] 0.156] 0.171] 0.156]
TYPE 2 RIT (reference = 1) SU-ARG (2) MgCl2-ARG (3) MgCl2-199 (4)
10.87 11.60 11.40 11.12
[10.69; [11.48; [11.63; [11.28;
11.05] 11.71] 11.86] 11.43]
10.98 11.72 11.61 11.19
0.00089 0.00183 0.00186 0.00117
[0.00029; [0.00131; [0.00080; [0.00061;
0.00148] 0.00236] 0.00293] 0,00174]
0.160 0.139 0.130 0.142
[0.143; [0.147; [0.146; [0.155;
0.178] 0.154] 0.161] 0.168]
TYPE3 RIT (reference = 1) SU-ARG (2) MgCl2-ARG (3) MgCl2-199 (4)
12.97 13.31 13.37 12.14
[12.83; [13.15; [13.19; [11.94;
13.12] 13.47] 13.56] 12.34]
12.85 13.33 13.54 12.20
0.00103 0.00133 0.00158 0.00180
[0.00053; [0.00079; [0.00088; [0.00086;
0.00152] 0.00188] 0.00229] 0.00274]
0.159 0.158 0.154 0.145
[0.146; 0.171] [0.147 0.169] [0.142; 0.166 [0.131; 0.159]
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
3579
Table 4 Results of the comparisons of parameter sets (h, i) between strains, performed with the F-testa Hypothesis — stability
Reference (1) = SUC-ARG (2) formulation Reference (1) = MgCl2-ARG (3) formulation Reference (1) =MgCl2-199 (4) formulation
Strain 1
Strain 2
Strain 3
F value
P value
F value
P value
F value
P value
43.1 43.6 30.0
1×10−10*** 9×10−11*** 1×10−8***
34.4 16.3 2.5
2×10−9*** 7×10−6*** 0.09 No different
44.8 46.8 3.9
6×10−11*** 3×10−11*** 0.03*
a * means a significant difference (0.01BPB0.05) between the parameter sets, and then a significant difference between the inactivation rate of the formulations. *** means a very high significant difference (PB0.0001) between the parameter sets, and then a very high significant difference between the inactivation rate of the formulations.
with U=(ln N0, h, i) for Eq. (3) or U =(h, i) for Eq. (4).
3.2. Model fitting All the data processing was performed on a Macintosh G3 working under Mac OS 8.1 with Mathematica® 3.0.1 (WOLFRAM Research, Inc., USA). The non-linear Eq. (3) was fitted to data by non-linear regression using the routine Nonlinear Fit of Mathematica. This iterative procedure of fit is based on a modified Levenberg–Marquardt algorithm [16]. It minimizes the sum of squares of the differences between calculated and observed values. The algorithm thus allows estimation of the parameter set corresponding to the lowest residual sum of squares (RSS) expressed as follows: RSS = [yi − yˆi ]2 that is n
RSS = % [ln N(t,T),i −(hti exp(iTi ))]2 i=1
where n is the number of experimental data, and ln N(t,T)i the ith measure of titer, ti the storage time of the ith measure and Ti the temperature of the ith measure. Initial values are required for parameters to begin the iterative procedure of estimation. In our case, initial values were estimated by using Eqs. (1) and (2), i.e. by fitting successively Equation I and II with the experimental data. When convergence of the estimation of the parameters with Eqs. (3) or (4) to reasonable values is reached (i.e. when the 2 merit function no longer decreases). The iterative process is stopped and the program gives the estimated parameters, with their asymptotic standard errors (obtained by linear approximation) and the asymptotic t-ratio calculated as the ratio of the parameter estimate to the approximate standard error. If the t-test is significant (P-value B 0.05), this means that the parameter estimate is significantly different from 0.
3.3. E6aluation of fit Two plots were used to assess the quality of fitting. First, the plot of standardized residuals (fit residuals scaled by their asymptotic standard errors) versus fitted values, is used to examine if residuals have a constant variance. Second, a quantile–quantile plot of standardized residuals versus a normal distribution which gives a direct check on the assumption of normality and goodness of fit and may also reveal suspicious points [17]. With non-linear regression, it is sometimes possible to converge to not consistent parameters, which are sensitive to the choice of the initial values. These values correspond in fact to a local minimum of the RSS function. This can be avoided in our program by running the iterative process several times, with various initial conditions randomly chosen in a realistic parameter space. If each run leads to the same parameter estimates, they can be considered as the optimal values, corresponding to the lowest RSS value, i.e. to a global minimum. The parameter correlation was appreciated with the asymptotic correlation matrix. If the correlation between two parameters is closed to 0.9, this means that the estimation of one of the parameters is almost entirely determined by the estimation of the other one [17]. The estimations are thus, less precise. Moreover, the diagonal terms of the asymptotic covariance matrix were used to estimate the marginal confidence parameter intervals: qi = q. i 9 t(h/2; n − p)(rii )1/2
(5)
where qi is the expected value of one of the parameters (ln N0, h or i); q. i is the estimated value of qi ; t(a/2; n − p) corresponds to the upper quantile of the Student’s distribution (h= 5%) with (n− p) degrees of freedom; n is the number of data and p the number of parameters; rii is the diagonal term of the asymptotic covariance matrix corresponding to qi. By linear approximation, a 95% confidence band can be estimated for a fixed T value, as the set of points for T fixed at a given value (ti, T, ln N(ti,T)) such as y0 = f(t0,T)9 (p× F(1 − h; p,n − p) × varf(t0,T))1/2
(6)
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
3580
with F(1 − h; p,n − p) the quantile of the Fisher’s distribution (h = 5%) var[ f(t0T)]=|ˆ 2r G0(J TJ) − 1G T0 RSSmin n− p (f(t0, T, U) (f(t0, T, U) (f(t0, T, U) , G0 = , lnN0 h i
| 2r =
G0 =[1,−exp(iT)t0,− h exp(iT)t0]
n
n
|ˆ 2r (J TJ) − 1, is equivalent to the covariance matrix where (f(ti, T, U) J is the Jacobian matrix = (qj q = q.
3.4. Comparison of data sets The statistical method used to compare various data sets is based on a comparison of two residual sums of squares [18]. 1. The RSS estimated with a full model when one parameter set is estimated for each data set (RSSf). 2. The RSS estimated with a partial model when a common parameter set is estimated for all data sets (RSSp). The statistical test consists in n− pf RSSp − RSSf Fratio(full/partial) = × (7) pf − pp RSSf
Fig. 4. Filled contour plot of the global model (calculation of ln N0 −ln N with Eq. (4)) where we reported the pairs of parameter estimated with every data set. Bold numbers correspond to various formulations, which are compared; their position on the graph depends on the (h,i) values; the nearest to the dark surface the number is the lowest the loss at 20°C.
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
3581
where n is the total number of data; pf the number of parameters in the full model; and pp the number of parameters in the partial model. The null hypothesis corresponds to a no significant difference between the full and partial models. The Fratio is then tested against the Fisher value Ftheo = F(1− h; pc −pp; n− pc). If (Fratio BFtheo), the test is non-significant and the parameter set is the same for all the data sets. It is noteworthy that the Fratio is only approximated because it does not have an exact Fisher distribution for a non-linear model [18]. The F-test was used in this study to compare the vaccine stability of a new formulation to a reference standard. In this case, parameters lnN0 are fixed according to their estimated values and only parameters h and i were compared. If the test is not significant, we conclude that the stability of the vaccine with a new formulation is not significantly different from the stability of the reference product. We also plotted the curves f(t,T,U) = constant with U =(h,i) for fixed values of t and T (t =19 days; T =20°C), in the plane (h,i). The couples of parameter estimations (h,i) calculated with each data set were reported in this plot. The parameter h varied between 0 and 0.002 and i between 0 and 0.20. On this graph the titer loss varied between 0 and 2 ln CCID50 (Fig. 1). The lowest loss was obtained when (h,i) were near 0, corresponding to the darkest area CCID50 (Fig. 1). The highest loss was obtained when (h,i) were near (0.002; 0.2), corresponding to the white area CCID50 (Fig. 1).
the expected values from the experiments) and the number of iterations was not too high (Table 1). The model was fitted to the data sets of Leal et al. [12], Table 3 shows the resulting parameter values. There was no problem to estimate the parameters and their 95% confidence limits. The F-test was performed to compare the stability of each new vaccine formulation to the reference vaccine (number 1). The results are shown in Table 4. The partial model (P\ 0.05) describes the inactivation of strain 2 with stabilizer 1 (reference) or 4 (vaccine with a new formulation). This means that the thermostability of these two vaccines is not significantly different (h= 0.05). For all the other comparisons, the F-test value was below 0.05, a significant difference was detected between the different formulations. We reported in Fig. 4 that all the pairs of parameters (h,i) were estimated with every data set formulation on the graph f(t,T,U) = constant, with U= (h,i) for the fixed values of t and T (t= 19 days; T= 20°C), in the plane (h,i). The closer the cross sign was to the dark area, the lower was the loss for a conservation 19 days at 20°C. We observed that the reference (number 1) was near the dark area for the three viral strains, meaning that the minimum losses are obtained through this formulation.
4. Results
k= a exp(b/T)
4.1. Model fit to the data The model (Eq. (3)) was fitted to Data sets 1–4. Results are summarized in Table 1. All the parameters were significantly different from 0 (the P-value of the t-test is lower than 0.05). With all the data set we had a good fit of the model. As an example, the fit of the model to the Data set 2 (strain 1) is shown in Fig. 2, for 5, 25, 37 and 45°C. The accuracy of fit for data set 2 (strain 1) is shown in Figs. 2 and 3a and b. Residuals have a constant variance between −2 and 2 (Fig. 3a), except for one point. On the quantile– quantile plot (Fig. 3b) the repartition of the points are along a straight-line indicating that residuals have a normal distribution. A low correlation between ln N0 and h or i and a high positive correlation between h and i is deduced from the approximate correlation matrix (Table 2). Despite this correlation the estimation of the parameters was always possible (results were consistent with
5. Discussion As a first approach we try to use the non-linear Arrhenius equation (8)
where k is the slope estimated with Eq. (3), a and b the parameters of the model and T the temperature expressed in Kelvin, however, after 30 iterations the estimation of the parameters was not reached. On the other hand our model provided an accurate estimation after only 8 iterations (Table 2). The stability of live viral vaccines cannot be treated like the stability of classical non-biological drugs. Drugs are pure chemical products, while biological products often are a complex and structured assembly of proteins, saccharides, amino acids. This certainly explains why the non-linear Arrhenius equation did not fit into our experimental data. Our study demonstrated a good fit of the global model in the 5–50°C temperature range. With our model, the parameter estimation was possible with small sized data sets (Table 2). However, the confidence limit estimates may be affected by many factors, such as the number of samples, sampling times and the choice of temperatures. Moreover, the sampling points have ideally to be equally spaced in the range studied.
3582
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
Finally, the variability of the measure can be important in biology and must be taken into account in the model; when several independent measures exist, all of the data must be used instead of the mean. Whatever the data set, we observed a high correlation between parameters h and i. This is indeed a structural correlation between the two parameters. It was attributed to the nature of the mathematical function that correlates them. As a consequence, a change in the estimate of h may be compensated by a simultaneous change in the estimate of i. We observed that the presence of this correlation did not affect the estimation of the parameters. The loss is minimal when the couple of parameters (h,i) have the smaller values (Fig. 3). With our experimental data low values of h were not especially associated with low values of i (Tables 1 and 2). Therefore, no interpretation may be done with only one parameter h or i. The parameter set (h,i) must be analyzed as a whole. The model is useful to compare the effect of different formulations. The F-test is a statistical test to discriminate the effect of the formulations on the stability. In this study we only compared two products, but the test can be applied to an unlimited number of products. The curves f(t,T,U) =constant with U= (h,i)provide a visual information on the stability of the different products and allow to rank formulation for their capacity to prevent inactivation. They are thus a helpful tool to sort the effect of a formulation on the viral inactivation. Our method should save time and provide rational in decision-making concerning choices among the development or production of vaccines. The advantage of our non-linear approach is in its use of viral titer, time and temperature to provide a direct estimation of the parameters. Our model could be a solution when a fit to the Arrhenius equation is not possible. It could be used to compare the stability of viruses suspended in different media. Moreover, this method integrates the biological variability of the titration, which is often important for the biological products. According to the thermal sensitivity of the specific virus, it is important to evaluate the fitting of the model before using it for comparisons. The model described in this paper provides a good description of the inactivation of virus stored in the liquid form at different temperatures. The comparison tests provide a convenient method of quickly comparing a new product at the developmental stage or products at different stages of the production to a reference.
Acknowledgements We thank Janine Viret, Colette Joyet and Laure Barjhoux for technical assistance and Marie Laure Delignette and Bernard Meignier for helpful discussions. We gratefully acknowledge the help of M. Luciani and B. Carpick for proofreading the manuscript.
References [1] Guardabasso V, Munson PJ, Rodbard D. EXPFIT: a program for simultaneous analysis of families of exponential decay curves. Comput Methods Progr Biomed 1988;27:55 – 63. [2] Peetermans J. Factors affecting the stability of viral vaccines. Dev Biol Stand 1996;87:87 – 101. [3] Hurst CJ, Goyke T. Survival of indigenous enteric viruses during storage of waste water sludge samples. Can J Microbiol 1986;32:645 – 8. [4] Biziagos E, Passagot J, Crance JM, et al. Long-term survival of hepatitis A virus and poliovirus type 1 in mineral water. Appl Environ Microbiol 1988;54:2705 – 10. [5] Lopes OS, Guimaraes SS, de Carvalho R. Studies on yellow fever vaccine. II Stability of the reconstituted product. J Biol Stand 1988;16:71 – 6. [6] Hurst CJ, Benton WH, McClellan KA. Thermal and water source effects upon the stability of enteroviruses in surface freshwaters. Can J Microbiol 1989;35:474 – 80. [7] Ishak R, Howard CR. The thermal stability of yellow fever vaccines. Memorias Do Institudo Oswaldo 1990;85:339 –45. [8] Depner K, Bauer T, Liess B. Thermal and pH stability of pestiviruses. Revue Scientifique et Technique 1992;11:885 –93. [9] Graham ML, Springthorpe VS, Sattar SA. Ex vivo protocol for testing virus survival on human skin: experiments with herpesvirus 2. Appl Environ Microbiol 1996;62:4252 – 5. [10] Kaplan C. The heat inactivation of Vaccinia Virus. J Genetic Microbiol 1958;18:58 – 63. [11] Mauler R, Gruschkau H. On stability of oral poliovirus vaccines. Dev Biol Stand 1978;41:267 – 70. [12] Leal ML, Lopes FJ, Carvalho MH, et al. Study on thermostabilizers for trivalent oral poliomyelitis vaccine. Memorias Do Institudo Oswaldo Cruz 1990;85:329 – 38. [13] King S-YP, Kung M-S. Statistical prediction of drug stability based on nonlinear parameter estimation. J Pharmaceut Sci 1984;73:657 – 62. [14] Ertel KD, Cartensen JT. Examination of a modified Arrhenius relationship for pharmaceutical stability prediction. Int J Pharmaceut 1990;61:9 – 14. [15] Husson-van Vliet J, Lemoine PE. Pre´ paration et e´ valuation d’une re´ fe´ rence interne de vaccin polio oral. Arch B Med Soc Hyg 1989;47:270 – 1. [16] Meyer RR, Roth PM. Modified damped least squares: an algorithm for nonlinear estimation. J Inst Math Appl 1972;9:218 – 33. [17] Bates D.M., Watts D.G. Nonlinear regression analysis and its applications, 1988. [18] Guardabasso V, Munson PJ, Rodbard D. EXPFIT: a program for simultaneous analysis of families of exponential decay curves. Comput Methods Progr Biomed 1988;27:55 – 63.
A mathematical model describing the thermal virus inactivation Catherine Noe¨l *, Sandrine Charles, Alain Franc¸on, Jean Pierre Flandrois De6elopment Department, A6entis Pasteur, 69280 Marcy l’Etoile, France Received 1 May 2000; received in revised form 4 December 2000; accepted 12 December 2000
Abstract A new mathematical model is proposed to describe the inactivation of viruses at different temperatures. This model takes into account the exponential decrease of the viral titer with time, the inactivation rate being an exponential function of the temperature. A one-step non-linear regression was used to fit oral poliovirus vaccine (OPV) experimental data. In one of the applications of the model, we illustrate the use of our model to compare the accelerated degradation test of OPV new formulations to standard OPV. Such a model is both simple and convenient to use. It should be a useful tool in optimizing formulations for live viral vaccines. © 2001 Elsevier Science Ltd. All rights reserved. Keywords: Kinetic; Mathematical model; Oral polymyelitis vaccine; Temperature; Thermal inactivation; Viral vaccines
1. Introduction Stability of biological products is a major requirement in the pharmaceutical industry [1], ensuring quality, safety of use and, for vaccines success of the immunization. Temperature is the major environmental factor that affects vaccine stability [2]. It is studied because of its practical importance as a limiting factor of the stability of live viral vaccine against polio (OPV), measles, mumps, rubella, varicella. The common method used to study the thermostability of live virus vaccines is to analyze the loss of titer (i.e. the decrease in virus number) with time, during storage at elevated temperatures for which inactivation will occur more rapidly. This is the so-called accelerated inactivation test or accelerated degradation test. The test usually consists of two steps (i) description of the exponential decrease of the titer with time [3–9], by computing the inactivation rate versus time after linearization; and (ii) the use of the Arrhenius equation which correlates the natural logarithm of the inactivation rate with the inverse of temperature [7,10 – 12].
* Corresponding author. Tel.: + 33-4-37373456; fax: + 33-437373029. E-mail address: [email protected] (C. Noe¨l).
This two-step approach is up to now not satisfactory; the error associated with the determination of viral titer are not taken into account, and therefore, the confidence interval parameters may be too wide to apply the model in a reliable predictive manner [13,14]. A more sophisticated way is to establish a mathematical relationship to directly describe the decrease of the viral titer as a function of both time and temperature during storage. We previously tried to build such a model based on the Arrhenius equation under its exponential form but it was impossible to fit the model to our data sets. The Arrhenius equation is based on a determinist mechanistic approach. On the other hand, we propose in this work an empirical approach with a formulation of the model based on simple mathematical functions, chosen for their simplicity and their suitability to fit our data. Our aim was to accurately describe the biological phenomenon with a robust and simple model, which can be used, for example, to compare new formulations of live viral vaccines. We, thus, developed this model on oral polioviruse vaccine and we applied it to compare the stability of different formulations. This paper focuses first on the non-linear fit of the mathematical function to experimental data (either from bioassays or literature), and second on the use of the model in comparison with virus inactivation in various formulations.
0264-410X/01/$ - see front matter © 2001 Elsevier Science Ltd. All rights reserved. PII: S0264-410X(01)00011-1
C. Noe¨l et al. / Vaccine 19 (2001) 3575–3582
3576
2. Materials and methods
2.1. Data We studied the inactivation of the three strains of Polio Sabin virus included in the oral polio vaccine (OPV), all of which are stored in liquid form. The data were either taken from the literature or obtained in Aventis Pasteur laboratories. Data set 1 was used to fit the model. It was extracted from Husson et al. [15]. The authors describe a trivalent
Table 1 Sampling schedule for Data set 2 and 3 to study the inactivation of polioviruses Data set
T (°C)
Sampling time
2
5 25 37 45 5 25 37
1, 3, 6, 9, 12, 18, 24, 30, 36 months 2, 4, 6, 8, 10, 12 weeks 2, 5, 7, 14, 21 days 2, 5, and 7 days 1, 3, 5, 6, 7 months 2, 4, 6, 8, 10, 12, 16 weeks 5, 7, 9, 10, 11, 12, 14, 16, 18, 19, 20, 21 days 2, 4, 6, 8, 10 days 1, 2, 3, 4, 5, 6 days 0.62, 0.83, 1.06, 1.25, 1.67 days
3
40 45 50
working reference preparation of OPV containing the three Sabin Strains. The global titer was measured after storage at 20°C for 1, 2, and 3 weeks; at 37°C for 1, 3 and 5 days; and at 45°C for 1, 2 and 3 days. The titer given in the table of that paper is the mean of two measures. An estimation of the titer at the beginning of the study time (t= 0) is given by the mean of 68 measures of the reference stored at − 70°C. Data sets 2 and 3 were used to fit the model. They were obtained in Aventis Pasteur laboratories. They provided the attenuated poliovirus strains used in this investigation: type 1: LSC 2ab/KP2 So + 2; type2: P712 ch 2ab/KP3 So +2; type3: Leon 12a1b/KP4 SO +4. In Data set 2, the titer was studied for OPV exposed to 5, 25, 37 or 45°C for 1–36 months (Table 1). Data set 3 concerns OPV exposed to six temperatures ranging for 5–50°C for 1–7 months (Table 1). Viral titer was defined in CCID50 per dose of vaccine. The titer of each strain was measured after neutralization of the two other strains using the appropriate antisera. Serial dilutions of virus suspension were performed using the culture medium, each dilution was inoculated in ten wells of a microplate and HepII cells were added in each well. The microplates were incubated at 32°C with CO2 and the cytopathogenic effect observed to determine the virus titer. Calculations were done using the Spearman–Karber method.
2.2. Comparison data set Data sets 4 used for comparison of various formulations were extracted from Leal et al. [12]. Different formulations of trivalent OPV were tested, each constituting a data set. The titer was analyzed after 4, 7, 14, 21, 28, 35 and 42 days at 20°C. 4, 7, 14, 21, 28 and 35 days at 26°C. 1, 2, 4, 7, 14 and 21 days at 37°C. The three virus strains were studied. Titers given in tables of the article represent the mean of four measures.
3. Modelling
3.1. Model building The inactivation of biological product exposed to heat is usually described by the following model [3–9]: Fig. 1. Filled contour plot of the global model Eq. (6) describing the effect of (h,i) on the titer loss (lnN0 − lnN), for a temperature of 20°C and a sampling time of 19 days: the titer losses are increasing from the black area to the white area. For example if h= 0.0005 and i = 0.05 the estimated titer loss is 0.02 log DICC50/ml and if h = 0.002 and i = 0.20 the estimated titer loss is 2.07 DICC50/ml when the vaccine is stored 19 days at 20°C. Moreover, all the points placed on the same line had an equal predicted loss.
Nt = N0 exp(− kt) U log Nt = log N0 − kt U ln Nt =ln N0 − kt
(1)
where Nt refers to the titer at time t, N0 to the initial titer, and k to the inactivation rate for a given temperature T. k has a biological interpretation as the inactivation rate (i.e. the loss of titer per unit of time).
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
3577
Fig. 2. Fit of the new global model Eq. (3) with Data set 2. Squares represent observed values, the continue-line the adjusted global model and the dash lines the 95% confidence limits calculated with Eq. (5).
Fig. 3. (a) Standardized residuals versus the calculated values of Data set 2; such a plot corroborates the good fit of the model; no obvious heteroscedasticity can be observed and the residuals seem to be randomly distributed. (b) Normal quantiles versus standardized residual quantiles of Data set 2; this plot shows the normality of the residual distribution and the absence of outliers.
To build the model, we studied the relation between k and the temperature T. We then propose a new relationship between k and T expressed as follows: k = h exp(iT)
(2)
where h correspond to the inactivation rate for temperature equal to 0, and i have no biological signification.
By substitution of Eq. (2) into Eq. (1), we finally get a complete model. ln N(t,T) = f(t, T, U) ln N(t,T) = ln N0 − h exp(iT)t
(3)
or ln N0 − ln N(t, T) = h exp(iT)t
(4)
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
3578
Table 2 Asymptotic results obtained by fitting model IV to Data sets 1, 2 and 3: parameter estimations, their 95% standard error, the t ratio and the P-value, and the correlation matrixa Fit data set
Virus titration
Parameter
Estimate
Standard Error
1 (n= 12)
Global
ln N0 h i
12.09 0.00800 0.127
ln N0 h i
2 (n = 81)
2 (n= 81)
2 (n=81)
3 (n= 146)
3 (n=134)
3 (n=149)
a
Strain: 1
Strain:2
Strain: 3
Strain: 1
Strain: 2
Strain: 3
t ratio
P value
Correlation matrix
Iteration number
0.26 0.00333 0.009
46.4 −2.4 13.9
B0.01 0.04 B0.01
Á 1 Ã Ã 0.58 Ã Ä−0.47
1 −0.98
 à à à 1Å
18
14.45 0.00080 0.149
0.10 0.00011 0.003
152.1 −7.4 47.3
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.61 Ã Ä−0.41
1 −0.93
 à à à 1Å
6
ln N0 h i
12.41 0.00065 0.150
0.10 0.00011 0.004
123.5 −5.9 37.7
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.61 Ã Ä−0.42
1 −0.94
 à à à 1Å
5
ln N0 h i
14.19 0.00085 0.141
0.10 0.00013 0.004
138.8 −6.8 39.2
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.60 Ã Ä−0.31
1 −0.89
 à à à 1Å
7
ln N0 h i
13.75 0.00083 0.163
0.11 0.00008 0.002
127.1 −10.3 69.1
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.31 Ã Ä−0.08
1 −0.96
 à à à 1Å
8
ln N0 h i
10.70 0.00118 0.154
0.14 0.00014 0.003
78.2 −8.7 56.2
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.35 Ã Ä−0.07
1 −0.94
 à à à 1Å
8
ln N0 h i
13.40 0.00107 0.153
0.15 0.00015 0.004
91.4 −7.0 43.0
B0.01 B0.01 B0.01
Á 1 Ã Ã 0.27 Ã Ä−0.02
1 −0.96
 à à à 1Å
7
The iteration number needed to estimate various parameters is also indicated.
Table 3 Asymptotic results obtained by fitting model IV to data sets 4: parameter estimations and their 95% standard error h Calculated
lnN0
i Calculated
Calculated
Observed
TYPE 1 RIT (reference = 1) SU-ARG (2) MgCl2-ARG (3) MgCl2-199 (4)
13.71 13.83 13.80 13.96
[13.83;13.95] [13.70; 13.95] [13.62; 13.97] [13.79; 14.12]
13.91 13.84 13.86 14.02
0,00079 0.00180 0.00133 0.00196
[0.00044; [0.00126; [0.00073; [0.00122;
0.00113] 0.00234] 0.00192] 0.00270]
0.166 0.149 0.159 0.146
[0.155; [0.140; [0.147; [0.135;
0.178] 0.156] 0.171] 0.156]
TYPE 2 RIT (reference = 1) SU-ARG (2) MgCl2-ARG (3) MgCl2-199 (4)
10.87 11.60 11.40 11.12
[10.69; [11.48; [11.63; [11.28;
11.05] 11.71] 11.86] 11.43]
10.98 11.72 11.61 11.19
0.00089 0.00183 0.00186 0.00117
[0.00029; [0.00131; [0.00080; [0.00061;
0.00148] 0.00236] 0.00293] 0,00174]
0.160 0.139 0.130 0.142
[0.143; [0.147; [0.146; [0.155;
0.178] 0.154] 0.161] 0.168]
TYPE3 RIT (reference = 1) SU-ARG (2) MgCl2-ARG (3) MgCl2-199 (4)
12.97 13.31 13.37 12.14
[12.83; [13.15; [13.19; [11.94;
13.12] 13.47] 13.56] 12.34]
12.85 13.33 13.54 12.20
0.00103 0.00133 0.00158 0.00180
[0.00053; [0.00079; [0.00088; [0.00086;
0.00152] 0.00188] 0.00229] 0.00274]
0.159 0.158 0.154 0.145
[0.146; 0.171] [0.147 0.169] [0.142; 0.166 [0.131; 0.159]
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
3579
Table 4 Results of the comparisons of parameter sets (h, i) between strains, performed with the F-testa Hypothesis — stability
Reference (1) = SUC-ARG (2) formulation Reference (1) = MgCl2-ARG (3) formulation Reference (1) =MgCl2-199 (4) formulation
Strain 1
Strain 2
Strain 3
F value
P value
F value
P value
F value
P value
43.1 43.6 30.0
1×10−10*** 9×10−11*** 1×10−8***
34.4 16.3 2.5
2×10−9*** 7×10−6*** 0.09 No different
44.8 46.8 3.9
6×10−11*** 3×10−11*** 0.03*
a * means a significant difference (0.01BPB0.05) between the parameter sets, and then a significant difference between the inactivation rate of the formulations. *** means a very high significant difference (PB0.0001) between the parameter sets, and then a very high significant difference between the inactivation rate of the formulations.
with U=(ln N0, h, i) for Eq. (3) or U =(h, i) for Eq. (4).
3.2. Model fitting All the data processing was performed on a Macintosh G3 working under Mac OS 8.1 with Mathematica® 3.0.1 (WOLFRAM Research, Inc., USA). The non-linear Eq. (3) was fitted to data by non-linear regression using the routine Nonlinear Fit of Mathematica. This iterative procedure of fit is based on a modified Levenberg–Marquardt algorithm [16]. It minimizes the sum of squares of the differences between calculated and observed values. The algorithm thus allows estimation of the parameter set corresponding to the lowest residual sum of squares (RSS) expressed as follows: RSS = [yi − yˆi ]2 that is n
RSS = % [ln N(t,T),i −(hti exp(iTi ))]2 i=1
where n is the number of experimental data, and ln N(t,T)i the ith measure of titer, ti the storage time of the ith measure and Ti the temperature of the ith measure. Initial values are required for parameters to begin the iterative procedure of estimation. In our case, initial values were estimated by using Eqs. (1) and (2), i.e. by fitting successively Equation I and II with the experimental data. When convergence of the estimation of the parameters with Eqs. (3) or (4) to reasonable values is reached (i.e. when the 2 merit function no longer decreases). The iterative process is stopped and the program gives the estimated parameters, with their asymptotic standard errors (obtained by linear approximation) and the asymptotic t-ratio calculated as the ratio of the parameter estimate to the approximate standard error. If the t-test is significant (P-value B 0.05), this means that the parameter estimate is significantly different from 0.
3.3. E6aluation of fit Two plots were used to assess the quality of fitting. First, the plot of standardized residuals (fit residuals scaled by their asymptotic standard errors) versus fitted values, is used to examine if residuals have a constant variance. Second, a quantile–quantile plot of standardized residuals versus a normal distribution which gives a direct check on the assumption of normality and goodness of fit and may also reveal suspicious points [17]. With non-linear regression, it is sometimes possible to converge to not consistent parameters, which are sensitive to the choice of the initial values. These values correspond in fact to a local minimum of the RSS function. This can be avoided in our program by running the iterative process several times, with various initial conditions randomly chosen in a realistic parameter space. If each run leads to the same parameter estimates, they can be considered as the optimal values, corresponding to the lowest RSS value, i.e. to a global minimum. The parameter correlation was appreciated with the asymptotic correlation matrix. If the correlation between two parameters is closed to 0.9, this means that the estimation of one of the parameters is almost entirely determined by the estimation of the other one [17]. The estimations are thus, less precise. Moreover, the diagonal terms of the asymptotic covariance matrix were used to estimate the marginal confidence parameter intervals: qi = q. i 9 t(h/2; n − p)(rii )1/2
(5)
where qi is the expected value of one of the parameters (ln N0, h or i); q. i is the estimated value of qi ; t(a/2; n − p) corresponds to the upper quantile of the Student’s distribution (h= 5%) with (n− p) degrees of freedom; n is the number of data and p the number of parameters; rii is the diagonal term of the asymptotic covariance matrix corresponding to qi. By linear approximation, a 95% confidence band can be estimated for a fixed T value, as the set of points for T fixed at a given value (ti, T, ln N(ti,T)) such as y0 = f(t0,T)9 (p× F(1 − h; p,n − p) × varf(t0,T))1/2
(6)
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
3580
with F(1 − h; p,n − p) the quantile of the Fisher’s distribution (h = 5%) var[ f(t0T)]=|ˆ 2r G0(J TJ) − 1G T0 RSSmin n− p (f(t0, T, U) (f(t0, T, U) (f(t0, T, U) , G0 = , lnN0 h i
| 2r =
G0 =[1,−exp(iT)t0,− h exp(iT)t0]
n
n
|ˆ 2r (J TJ) − 1, is equivalent to the covariance matrix where (f(ti, T, U) J is the Jacobian matrix = (qj q = q.
3.4. Comparison of data sets The statistical method used to compare various data sets is based on a comparison of two residual sums of squares [18]. 1. The RSS estimated with a full model when one parameter set is estimated for each data set (RSSf). 2. The RSS estimated with a partial model when a common parameter set is estimated for all data sets (RSSp). The statistical test consists in n− pf RSSp − RSSf Fratio(full/partial) = × (7) pf − pp RSSf
Fig. 4. Filled contour plot of the global model (calculation of ln N0 −ln N with Eq. (4)) where we reported the pairs of parameter estimated with every data set. Bold numbers correspond to various formulations, which are compared; their position on the graph depends on the (h,i) values; the nearest to the dark surface the number is the lowest the loss at 20°C.
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
3581
where n is the total number of data; pf the number of parameters in the full model; and pp the number of parameters in the partial model. The null hypothesis corresponds to a no significant difference between the full and partial models. The Fratio is then tested against the Fisher value Ftheo = F(1− h; pc −pp; n− pc). If (Fratio BFtheo), the test is non-significant and the parameter set is the same for all the data sets. It is noteworthy that the Fratio is only approximated because it does not have an exact Fisher distribution for a non-linear model [18]. The F-test was used in this study to compare the vaccine stability of a new formulation to a reference standard. In this case, parameters lnN0 are fixed according to their estimated values and only parameters h and i were compared. If the test is not significant, we conclude that the stability of the vaccine with a new formulation is not significantly different from the stability of the reference product. We also plotted the curves f(t,T,U) = constant with U =(h,i) for fixed values of t and T (t =19 days; T =20°C), in the plane (h,i). The couples of parameter estimations (h,i) calculated with each data set were reported in this plot. The parameter h varied between 0 and 0.002 and i between 0 and 0.20. On this graph the titer loss varied between 0 and 2 ln CCID50 (Fig. 1). The lowest loss was obtained when (h,i) were near 0, corresponding to the darkest area CCID50 (Fig. 1). The highest loss was obtained when (h,i) were near (0.002; 0.2), corresponding to the white area CCID50 (Fig. 1).
the expected values from the experiments) and the number of iterations was not too high (Table 1). The model was fitted to the data sets of Leal et al. [12], Table 3 shows the resulting parameter values. There was no problem to estimate the parameters and their 95% confidence limits. The F-test was performed to compare the stability of each new vaccine formulation to the reference vaccine (number 1). The results are shown in Table 4. The partial model (P\ 0.05) describes the inactivation of strain 2 with stabilizer 1 (reference) or 4 (vaccine with a new formulation). This means that the thermostability of these two vaccines is not significantly different (h= 0.05). For all the other comparisons, the F-test value was below 0.05, a significant difference was detected between the different formulations. We reported in Fig. 4 that all the pairs of parameters (h,i) were estimated with every data set formulation on the graph f(t,T,U) = constant, with U= (h,i) for the fixed values of t and T (t= 19 days; T= 20°C), in the plane (h,i). The closer the cross sign was to the dark area, the lower was the loss for a conservation 19 days at 20°C. We observed that the reference (number 1) was near the dark area for the three viral strains, meaning that the minimum losses are obtained through this formulation.
4. Results
k= a exp(b/T)
4.1. Model fit to the data The model (Eq. (3)) was fitted to Data sets 1–4. Results are summarized in Table 1. All the parameters were significantly different from 0 (the P-value of the t-test is lower than 0.05). With all the data set we had a good fit of the model. As an example, the fit of the model to the Data set 2 (strain 1) is shown in Fig. 2, for 5, 25, 37 and 45°C. The accuracy of fit for data set 2 (strain 1) is shown in Figs. 2 and 3a and b. Residuals have a constant variance between −2 and 2 (Fig. 3a), except for one point. On the quantile– quantile plot (Fig. 3b) the repartition of the points are along a straight-line indicating that residuals have a normal distribution. A low correlation between ln N0 and h or i and a high positive correlation between h and i is deduced from the approximate correlation matrix (Table 2). Despite this correlation the estimation of the parameters was always possible (results were consistent with
5. Discussion As a first approach we try to use the non-linear Arrhenius equation (8)
where k is the slope estimated with Eq. (3), a and b the parameters of the model and T the temperature expressed in Kelvin, however, after 30 iterations the estimation of the parameters was not reached. On the other hand our model provided an accurate estimation after only 8 iterations (Table 2). The stability of live viral vaccines cannot be treated like the stability of classical non-biological drugs. Drugs are pure chemical products, while biological products often are a complex and structured assembly of proteins, saccharides, amino acids. This certainly explains why the non-linear Arrhenius equation did not fit into our experimental data. Our study demonstrated a good fit of the global model in the 5–50°C temperature range. With our model, the parameter estimation was possible with small sized data sets (Table 2). However, the confidence limit estimates may be affected by many factors, such as the number of samples, sampling times and the choice of temperatures. Moreover, the sampling points have ideally to be equally spaced in the range studied.
3582
C. Noe¨ l et al. / Vaccine 19 (2001) 3575–3582
Finally, the variability of the measure can be important in biology and must be taken into account in the model; when several independent measures exist, all of the data must be used instead of the mean. Whatever the data set, we observed a high correlation between parameters h and i. This is indeed a structural correlation between the two parameters. It was attributed to the nature of the mathematical function that correlates them. As a consequence, a change in the estimate of h may be compensated by a simultaneous change in the estimate of i. We observed that the presence of this correlation did not affect the estimation of the parameters. The loss is minimal when the couple of parameters (h,i) have the smaller values (Fig. 3). With our experimental data low values of h were not especially associated with low values of i (Tables 1 and 2). Therefore, no interpretation may be done with only one parameter h or i. The parameter set (h,i) must be analyzed as a whole. The model is useful to compare the effect of different formulations. The F-test is a statistical test to discriminate the effect of the formulations on the stability. In this study we only compared two products, but the test can be applied to an unlimited number of products. The curves f(t,T,U) =constant with U= (h,i)provide a visual information on the stability of the different products and allow to rank formulation for their capacity to prevent inactivation. They are thus a helpful tool to sort the effect of a formulation on the viral inactivation. Our method should save time and provide rational in decision-making concerning choices among the development or production of vaccines. The advantage of our non-linear approach is in its use of viral titer, time and temperature to provide a direct estimation of the parameters. Our model could be a solution when a fit to the Arrhenius equation is not possible. It could be used to compare the stability of viruses suspended in different media. Moreover, this method integrates the biological variability of the titration, which is often important for the biological products. According to the thermal sensitivity of the specific virus, it is important to evaluate the fitting of the model before using it for comparisons. The model described in this paper provides a good description of the inactivation of virus stored in the liquid form at different temperatures. The comparison tests provide a convenient method of quickly comparing a new product at the developmental stage or products at different stages of the production to a reference.
Acknowledgements We thank Janine Viret, Colette Joyet and Laure Barjhoux for technical assistance and Marie Laure Delignette and Bernard Meignier for helpful discussions. We gratefully acknowledge the help of M. Luciani and B. Carpick for proofreading the manuscript.
References [1] Guardabasso V, Munson PJ, Rodbard D. EXPFIT: a program for simultaneous analysis of families of exponential decay curves. Comput Methods Progr Biomed 1988;27:55 – 63. [2] Peetermans J. Factors affecting the stability of viral vaccines. Dev Biol Stand 1996;87:87 – 101. [3] Hurst CJ, Goyke T. Survival of indigenous enteric viruses during storage of waste water sludge samples. Can J Microbiol 1986;32:645 – 8. [4] Biziagos E, Passagot J, Crance JM, et al. Long-term survival of hepatitis A virus and poliovirus type 1 in mineral water. Appl Environ Microbiol 1988;54:2705 – 10. [5] Lopes OS, Guimaraes SS, de Carvalho R. Studies on yellow fever vaccine. II Stability of the reconstituted product. J Biol Stand 1988;16:71 – 6. [6] Hurst CJ, Benton WH, McClellan KA. Thermal and water source effects upon the stability of enteroviruses in surface freshwaters. Can J Microbiol 1989;35:474 – 80. [7] Ishak R, Howard CR. The thermal stability of yellow fever vaccines. Memorias Do Institudo Oswaldo 1990;85:339 –45. [8] Depner K, Bauer T, Liess B. Thermal and pH stability of pestiviruses. Revue Scientifique et Technique 1992;11:885 –93. [9] Graham ML, Springthorpe VS, Sattar SA. Ex vivo protocol for testing virus survival on human skin: experiments with herpesvirus 2. Appl Environ Microbiol 1996;62:4252 – 5. [10] Kaplan C. The heat inactivation of Vaccinia Virus. J Genetic Microbiol 1958;18:58 – 63. [11] Mauler R, Gruschkau H. On stability of oral poliovirus vaccines. Dev Biol Stand 1978;41:267 – 70. [12] Leal ML, Lopes FJ, Carvalho MH, et al. Study on thermostabilizers for trivalent oral poliomyelitis vaccine. Memorias Do Institudo Oswaldo Cruz 1990;85:329 – 38. [13] King S-YP, Kung M-S. Statistical prediction of drug stability based on nonlinear parameter estimation. J Pharmaceut Sci 1984;73:657 – 62. [14] Ertel KD, Cartensen JT. Examination of a modified Arrhenius relationship for pharmaceutical stability prediction. Int J Pharmaceut 1990;61:9 – 14. [15] Husson-van Vliet J, Lemoine PE. Pre´ paration et e´ valuation d’une re´ fe´ rence interne de vaccin polio oral. Arch B Med Soc Hyg 1989;47:270 – 1. [16] Meyer RR, Roth PM. Modified damped least squares: an algorithm for nonlinear estimation. J Inst Math Appl 1972;9:218 – 33. [17] Bates D.M., Watts D.G. Nonlinear regression analysis and its applications, 1988. [18] Guardabasso V, Munson PJ, Rodbard D. EXPFIT: a program for simultaneous analysis of families of exponential decay curves. Comput Methods Progr Biomed 1988;27:55 – 63.