Abstract
Standardized statistical and graphical methods for analysis of limiting dilution assays are highly desirable to enable investigators to compare and interpret results and conclusions with greater accuracy and precision. According to these requirements, we present in this work a powerful statistical slope test that estimates the fit of the singlehit Poisson model to limiting dilution experiments. This method is readily amenable to a graphical representation. This slope test is obtained by modeling limiting dilution data according to a linear loglog regression model, which is a generalized linear model specially designed for modeling binary data. The result of the statistical slope test can then be graphed to visualize whether the data are compatible or not with the singlehit Poisson model. We demonstrate this statistical test and its graphical representation by using two examples: a real limiting dilution experiment evaluating the growth frequency of IL2responsive tumorinfiltrating T cells in a malignant lymph node involved by a B cell nonHodgkin’s lymphoma, and a simulation of a limiting dilution assay corresponding to a theoretical nonsinglehit Poisson model, suppressor twotarget Poisson model.
Limiting dilution analysis has gained widespread acceptance as a tool for quantifying the frequency of cells in the immune system that possess various functional activities (1). To analyze the data, it is usually assumed that only one limiting cell of only one cell subset is necessary and sufficient for generating a positive response (the singlehit Poisson model (SHPM)^{2}. However, the immune system is a complex network of cellular and humoral (cytokines) interactions. It is thus not unreasonable to consider that experimental data do not adhere to allornone functionality and will yield nonlinear titration curves (1, 2). In the literature, deviation from singlehit kinetics is usually estimated by a standard χ^{2} test. Previously, we demonstrated that the standard χ^{2} test was insufficient for estimating the goodnessoffit to the SHPM due to its intrinsic lack of power (3). Subsequently, we presented evidence that modeling limiting dilution assays by a loglog generalized linear model results in a statistical slope test that is better able to account for deviation from linearity than the standard χ^{2} test (4). In the present paper we show that, in addition to its ability to accurately detect departure of the experimental data from the SHPM, this slope test is readily amenable to a graphical representation.
In recent years, the continuing development of more complex models giving mathematical representation of biological data, in parallel with the tremendous improvements in computing power, have led to the emergence of even more complex statistical methodologies. In this context, generalized linear theorybased statistical tests represent a major advance in evaluating the fit of the SHPM to limiting dilution analysis. However, new methods are rarely used in practice, and a major challenge for research statisticians is to explain to immunologists the benefit that they can draw from these new methods. In line with this problem, it must be emphasized that graphing can make the result of a statistical investigation clearer to the reader and more attractive than if the result has been presented merely in a table. In this way, we present an understandable graphical illustration of a generalized linear modelbased statistical slope test with the aim of encouraging investigators to use a new statistical method for their limiting dilution data.
The generalized linear modelbased statistical procedure, accompanied by its graphical representation, is illustrated in this paper with two examples. The first set of data is derived from a real experiment assessing the growth frequency of IL2responsive tumorinfiltrating lymphocyte T (TILT) in a malignant lymph node involved by a lowgrade B cell nonHodgkin’s lymphoma. The second set of data is artificially generated and corresponds to a theoretical nonSHPM, suppressor twotarget Poisson model (STTPM). Specifically, in this second example, the total cell population contains a suppressor cell subset exerting a suppressor activity on the growth of a second cell subset (3). A possible physiologic illustration of this situation is given by the functional relationship which has been described between CD4^{+} virgin and memory T cells, where memory cells produce IL10, which in return inhibits the growth of virgin cells (2, 5).
Materials and Methods
Abbreviations used in the text, f, k, i, n, x, r, m, and μ, are defined as follows: f is the estimate of the cell frequency according to the SHPM hypothesis; k is the number of groups of replicate wells, each group labeled i, with i = 1, 2, 3, … , k. At each i, n_{i} is the number of replicate wells, x_{i} is the number of cells plated in each replicate well, r_{i} is the number of observed negative wells, m_{i} is the observed fraction of negative wells, m_{i} = r_{i}/n_{i}, μ_{i} is the theoretical fraction of negative wells.
Generalized linear model fitting the SHPM
The first step consists of modeling limiting dilution data according to a generalized linear loglog model fitting the SHPM and checking this model by an appropriate slope test (4). According to the zero term of the Poisson distribution, the theoretical fraction μ_{i} of negative wells is given by μ_{i} = exp(−fx_{i}).
This equation can be written according to a loglog linear model as log(−log(μ_{i})) = log(f) + log(x_{i}). This equation is now in the form Y_{i} = α + βX_{i} with Y_{i} = log(−log(μ_{i})), α = log(f), β = 1, and X_{i} = log(x_{i}).
A test of deviation from the SHPM (model checking) is a test of whether the estimate of the slope β is significantly different from 1. The appropriate statistical test z uses the normal deviates under the null hypothesis that β = 1. The hypothesis that the slope β is compatible with 1 at the 95% confidence level can also be expressed as This statement can be rewritten as β − 1.96 SE (β) < 1 < β + 1.96 SE (β)showing that the value, 1, is included in the 95% confidence interval of the slope β if the SHPM hypothesis holds.
Maximum likelihood estimates of the parameters α and β were obtained with the Fisher’s method of scoring according to an iteratively reweighted leastsquares (IRLS) procedure written for binary data (4, 6, 7). Maximum likelihood estimates satisfy the equation X^{T}WXβ = X^{T}WZ with X, the matrix with elements X_{i} = log (x_{i}); W, the diagonal matrix of weights with elements where V_{i} is the variance of μ_{i} and thus Z is the matrix of dependent variables with elements and thus β is the matrix of the coefficients α and β. The revised estimateβ^{(m)} solved iteratively by IRLS is β^{(m)} = (X^{T}WX)^{−1}X^{T}WZ where m is the m^{th} iteration and β^{(m)} is the matrix of the coefficients α^{(m)}and β^{(m)}. The elements α and β of the matrix β^{(m)} are given by The estimates of the variances for α and β, var(α) and var(β), are given bythe diagonal elements of the variancecovariance matrix (X^{T}WX)^{−1} evaluated with α^{(m)} β^{(m)}. The variance of β is given by The 95% confidence interval for β was calculated as The main steps of the IRLS algorithm used for calculations of α, β, and var(β) are given in the Appendix.
Computation of the cell subset frequency according to the SHPM
The second step consists in computation of the cell subset frequency according to the SHPM by maximum likelihood estimation. Let f be the estimate of the cell frequency; the maximum likelihood of f is the value of f that maximizes where log(L) is the natural logarithm of the likelihood function L (8), and P_{i} is given by P_{i} = exp(−fx_{i}) according to the SHPM.
The variance of f was calculated as the negative reciprocal of the second derivative of log(L) (8). The 95% confidence interval for f was calculated as 95% CI(f ) = f ± 1.96 SE (f ).
Derivation of the STPPM
For simulation experiments, it was assumed that the total CD4^{+} T cell population contains suppressor cells with frequency φ, demonstrating a growth suppressor activity, for instance, by the production of IL10, on another population of CD4^{+} T cells capable to growth with frequency γ (3). According to the Poisson law, the fraction of wells not containing suppressor cells, but containing at least one proliferating Tlymphocyte precursor cell, is exp(−φx_{i})(1 − exp(−γx_{i})). Thus, the theoretical fraction μ_{i} of negative wells is μ_{i} = 1 − (exp(−φx_{i})(1 − exp(−γx_{i}))). Calculations were made with φ = 0.08 and γ = 0.4. The theoretical number of negative wells was equal to 100 μ_{i} (100 replicate wells for each cell dose; n_{i} = 100), and the experimental number of negative wells, r_{i}, was taken at the nearest integer.
Results and Discussion
Table I⇓ presents the experimental and simulated limiting dilution data. Table II⇓ indicates the results of statistical analysis of the data and precursor frequency. In the real experiment conducted with TILT in one case of B cell lymphoma (9), the slope test z derived from the generalized linear loglog model did not reject the SHPM, and it was concluded that the SHPM adequately fits the observed data. In contrast, in the simulated experiment, conducted according to the STTPM, the statistics z clearly rejected the hypothesis that the data conformed to singlehit kinetics. This reflects the fact that the slope β was not compatible with 1 (95% confidence interval for the slope β ranged from 0.598 to 0.936).
Fig. 1⇓A and Fig. 2⇓ are graphical representations of the slope test z for both examples. This graphical representation is designed to allow visualization of whether the slope β of the loglog regression line fitted to the experimental data is compatible or not with 1 (SHPM adequacy or inadequacy). The graph is constructed in the following way. The loglog transform of the proportion of negative wells is plotted against the logarithm of the number of T cells per well. First, the observed data points (X_{i}, log(−log (m_{i})) are plotted. Then the loglog regression line fitted to the experimental data (regression line with equation Y_{i} = α + βX_{i}; Fig. 1⇓A and Fig. 2⇓, blue plain line) is plotted. Next are plotted two regression lines with the same intercept α as the loglog regression line, but with slope β equal to either the upper or the lower value of the 95% confidence interval for β (boundary regression lines, Y_{i} = α + β_{upper}X_{i}, Y_{i} = α + β_{lower}X_{i}; Fig. 1⇓A and Fig. 2⇓, blue dotted lines). Finally, a theoretical, socalled SHPM regression line, with the same intercept α, but with slope β equal to 1, is superimposed on the graph (Y_{i} = α + X_{i}; Fig. 1⇓A and Fig. 2⇓, red line). If the theoretical SHPM regression line with slope β equal to 1 does not lie between the two boundary regression lines for β_{upper} and β_{lower}, then the slope β of the loglog regression line fitted to the experimental data is not compatible with 1, and thus the SHPM does not fit the data. Reciprocally, if the SHPM regression line with slope β equal to 1 lies inside the limits of the two boundary regression lines for β_{upper} and β_{lower}, then the slope β of the loglog regression line fitted to the experimental data is compatible with 1, and thus the SHPM fits the data. To construct the graphs presented in Fig. 1⇓A and Fig. 2⇓, Table III⇓ gives the numerical values of the pairs (X_{i}, log(−log(m_{i})) for plotting observed data points, and Tables IV⇓ and V⇓ show determination of the pairs (X_{i}, Y_{i}) for drawing the four regression lines.
Fig. 1⇑A represents the graphical evaluation of the slope test z in the experiment evaluating the growth frequency of IL2responsive TILT in lymphoma. The equation of the loglog regression line fitted to the experimental data is Y_{i} = −3.081 + 0.909X_{i}, and the boundary equations of the two regression lines corresponding to the lower and upper limits of the slope β are Y_{i} = −3.081 +0.763X_{i} (lower value) and Y_{i} = −3.081 + 1.055X_{i} (upper value). Obviously, the regression line with β equal to 1, Y_{i} = −3.081 + X_{i}, lies inside the boundary regression lines for β = 0.763 and β = 1.055, thereby indicating that the experimental data are compatible with the SHPM. The graphical display of the slope test z in the STTPM is shown in Fig. 2⇑. The equation of the loglog regression line fitted to the experimental data is Y_{i} = −1.095 + 0.767X_{i}, and the boundary equations of the two regression lines for the lower and upper limits of the slope β are Y_{i} = −1.095 + 0.598X_{i} (lower value) and Y_{i} = −1.095 + 0.936X_{i} (upper value). The regression line with β equal to 1, Y_{i} = −1.095 + X_{i}, lies outside the boundary regression lines for β = 0.598 and β = 0.936, indicating that the data do not conform to the SHPM, thus precluding the computation of the cell frequency f according to singlehit kinetics.
Because the limiting dilution assay with TILT is compatible with the SHPM, results can be presented by the usual graphical display of limiting dilution data (semilog plot, Fig. 1⇑B). The fraction of negative wells on a logarithmic scale is plotted against the number of T cells per well. The prediction equation of the regression line is −log(μ_{i}) = 0.0339x_{i}, where 0.0339 is the cell frequency f estimated by the maximum likelihood method according to the SHPM hypothesis. Lower and upper dotted lines in Fig. 1⇑B are plotted by using lower and upper values of the 95% confidence interval of f. The boundary equations are −log(μ_{i}) = 0.0289x_{i} (lower), and −log(μ_{i}) = 0.0391x_{i} (upper).
Based on these findings, we propose that immunologists adopt a standard twograph method for presenting the results of limiting dilution data. The first graph should present the loglog transform of the results with the graphical representation of the slope test z designed to visualize whether the data are consistent or not with the SHPM. Four regression lines should be presented: the loglog regression line fitted to the experimental data according to the generalized linear loglog model, the two boundary regression lines for β_{upper} and β_{lower}, and the theoretical regression line with β equal to 1 (SHPM regression line). The graph should be accompanied by the equation of the loglog regression line fitted to the experimental data, the SE of β, the 95% confidence interval of β, and the numerical values of z and p(z) (Fig. 1⇑A and Fig. 2⇑).
If the SHPM hypothesis holds, the second graph should be the usual graphical representation of limiting dilution data (semilog plot), where the fraction of negative wells is plotted on a logarithmic scale against the number of cells per well. The prediction equation of the regression line is −log(μ_{i}) = fx_{i}, with the frequency f estimated by the maximum likelihood method according to the SHPM hypothesis. On this graph should also be plotted the two regression lines based on the extreme values of the confidenceinterval for f. The graph should be accompanied by the numerical value of f (or its inverse), its 95% confidence interval, and the equation of the regression line (Fig. 1⇑B).
In conclusion, this proposition broadly resembles a twograph standard presentation of data adopted a number of years ago by immunologists dealing with radioimmunoassays. The first graph eliminates nonspecific binding by demonstrating that the concentration of bound ligand (yaxis) increases nearly linearly with increasing free ligand (xaxis), then asymptotically approaches a plateau; this graph serves as a validation for transforming the data according to a second graph (the socalled Scatchard plot; bound/free vs bound) used to determine the number of binding sites and the ligand binding affinity. In the same manner, data from limiting dilution assays could adequately be summarized and validated according to the SHPM by routine representation using the two graphs presented in this report.
Appendix 1
Calculations of α, β, and var(β) according to the IRLS algorithm
A convenient feature of this algorithm is that it suggests a simple starting procedure to get the iteration under way. This consists of using the data themselves as the initial values of the dependent variable, Z_{0}, and the weight, W_{0}, and from these values are then calculated the first values of the parameter estimates, α_{0} and β_{0}. These first values for α and β are in turn used to calculate the first value of the theoretical fraction of negative wells, μ_{1}, which in turn is used to obtain the revised values of weight and dependent variable, Z_{1} and W_{1}, which in turn are used to calculate the revised values of the parameter estimates, α_{1} and β_{1}, and so on. The iteration is continued until the absolute change in the values of α_{m} and β_{m} in two successive cycles of the iterative process is <0.0001. m is the mth cycle of the iteration. The estimates of α and β, thus obtained, are denoted by α̂ and β̂, respectively. The corresponding fitted values of weights, W, are denoted by Ŵ.
First step: initiation of the algorithm and calculation of the starting values, W_{0}, Z_{0}, α_{0}, and β_{0}
For each group of replicate wells, i, 1) compute the initial value of the weight, W_{i0} 2) Compute the initial value of the dependent variable, Z_{i0} 3) Compute the first values of the parameter estimates, α_{0} and β_{0}, according to equations 4 and 5, using W_{i0} and Z_{i0}.
Second step, the iterative process for calculating the revised values of the parameter estimates, α_{m} and β_{m}
1) Compute the estimates of the theoretical fraction of negative wells, μ_{im}: μ_{im} = exp(−exp(α_{m − 1} + β_{m − 1}X_{i})). 2) Compute W_{im} according to equation 2 using μ_{im}. 3) Compute Z_{im} according to equation 3 using μ_{im} and α_{m − 1}, β_{m − 1}. 4) Compute the revised values of the parameter estimates, α_{m} and β_{m}, according to equations 4 and 5, using W_{im} and Z_{im}.
Repeat steps one to four until the absolute change in the values of α_{m} and β_{m} in two successive cycles of the iterative process is sufficiently small (<0.0001). The estimates of α and β are denoted by α̂ and β̂. The corresponding fitted values of weights are denoted by Ŵ_{i}.
Final step
1) Compute the variance of β̂ and var(β̂), according to equation (6), using the values Ŵ_{i}. 2) Compute the SE of β̂ 3) Compute the value, z, of the slope test according to equation 1 using β and
4) Compute the 95% confidence interval for β̂ according to equation (7) using β̂ and SE(β̂).
Acknowledgments
We thank J. L. Martiel (Unité Mixte de Recherche, Centre National de la Recherche Scientifique 5525), J. J. Lawrence (Institut National de la Santé et de la Recherche Médicale Unité 309), and E. Deslandres for helpful discussions throughout the development of this work.
Footnotes

↵1 Address correspondence and reprint requests to Dr. Thierry Bonnefoix, Groupe de Recherche sur les Lymphomes, Institut Albert Bonniot, RondPoint de la Chantourne, 38706 La Tronche, France. Email address: Thierry.Bonnefoix{at}ujfgrenoble.fr

↵2 Abbreviations used in this paper: SHPM, singlehit Poisson model; TILT, tumorinfiltrating lymphocyte T; STTPM, suppressor twotarget Poisson model; IRLS, iteratively reweighted leastsquares.
 Received February 8, 2001.
 Accepted September 18, 2001.
 Copyright © 2001 by The American Association of Immunologists