Abstract
Although there is currently no doubt that regulatory lymphocytes represent a master player in the immune system, a major unresolved problem is the accurate quantitation of these cells among unfractionated cell populations. This difficulty mainly arises because there are no specific immunophenotypic markers that can reliably discriminate between effector and regulatory lymphocytes. To face this problem, we have developed computational models of limiting dilution analyses addressing the question of the accurate estimation of the frequencies of effector and regulatory cells functionally engaged in an immune response. A set of generic equations were provided to form a framework for modeling limiting dilution data, enabling discrimination between qualitatively different models of suppression. These models include either one or two subpopulations of regulatory cells, featured by either low or potent regulatory activity. The potential of this modeling approach was illustrated by the accurate determination of the frequencies of effector and regulatory T lymphocytes in one real limiting dilution experiment of CD4^{+}CD25^{+} T lymphocytes performed in the context of an allogeneic response in the human system. The crucial advantage of the limiting dilution method over the “static, phenotypebased” method is the dynamic evaluation of effector and regulatory T cell biology through their actual functional activity.
In spite of the experimental evidence of the existence of regulatory lymphoid cells, the exact number of regulatory cell subpopulations, their differentiation pathways, and their functional relationships remain intense matters of debate (1, 2, 3, 4, 5, 6, 7, 8, 9). Studies designed to find molecular markers that are specific for regulatory lymphocytes are in progress in both mice and human systems, and the list of candidate molecules is growing. At present, the markers that best identify regulatory cells include CD25, CTLA4 (CD152), glucocorticoidinduced TNF receptor familyrelated gene (GITR^{3} or TNFRSF18), CD103, and, recently, neuropilin1 (10). Although it is possible to use these molecules to enrich for T cells with regulatory activity, none of them has been found that fulfills the criteria of exclusivity to regulatory lymphocytes and stability. In contrast, the expression of the transcription factor Foxp3 has been claimed to be a unique and essential attribute of naturally occurring, or activationinduced, CD4^{+}CD25^{+} regulatory T lymphocytes (11, 12, 13). However, this statement must be softened by several concerns suggesting that, although Foxp3 seems necessary for the generation of regulatory lymphocytes (14, 15, 16), it may not be sufficient regarding the suppressor function of these cells. First, it has not been unambiguously determined whether all Foxp3^{+}, scurfinexpressing T cells are actually capable of exerting potent regulatory activity, and a functional heterogeneity within this population cannot be excluded. Studies performed at the clonal level of Foxp3^{+} T cells should be decisive to resolve the question of the exact relationship between Foxp3 gene expression and regulatory activity. Second, we also have to consider the possibility that the suppressor activity of regulatory cells may heavily depend on both the nature and strength of their stimulation (17, 18, 19). Recently, a number of studies have demonstrated that external, costimulatory signals can inactivate regulatory T cells: GITR ligand, expressed by APC, and 41BB ligand (CD137 ligand) can abrogate the suppressive activity of CD4^{+}CD25^{+} GITR^{+} (20, 21, 22, 23) or CD4^{+}CD25^{+} 41BB^{+} regulatory cells (24). Third, Foxp3negative clonal cells derived from type1 (IL10 producing) regulatory T cells have been described (6). Recently, it has been unambiguously demonstrated that IL10secreting regulatory T cells, generated through in vitro or in vivo strategies, do no express the transcription factor Foxp3 (9, 25). Fourth, the cytokine IL6, produced by APC, can render pathogenspecific T cells refractory to the suppressive activity of CD4+CD25+ regulatory cells (26, 27). Fifth, at the opposite, other factors may be required for optimal regulatory function of these cells, like the STAT1 (28). Finally, even Th1 and Th2 cells, by their production of IL10, can clearly function as “regulatory” cells (9). As a direct consequence, a critical question remains: among unfractionated populations of lymphoid cells, what is the fraction of regulatory cells actually engaged in an immune response?
Taken together, these findings strongly encourage the idea that the biology of regulatory cells must be evaluated through their actual functional activity in the global context of their interaction with other cell subsets, rather than based on their phenotypic, surface or intracellular (Foxp3), markers. To face this situation, the limiting dilution method should be a tool of choice for accurately quantitating regulatory T cells regardless of their phenotypic markers and actually engaged in an immune response. Limiting dilution analysis has gained widespread acceptance as a standard method for quantifying accurately the frequency of cells that possess various functional activities and should be an invaluable tool to model the dynamics of functional interaction between cellular subpopulations (29, 30). Recently, we elaborated a regulatory cellbased mathematical model, describing the interaction between regulatory and proliferating T cells (31). In the present study, we have developed two additional models, which extend this previous model of T cellmediated suppression. These probabilistic models were applied to a real limiting dilution (LD) data set originating from a report on the suppressive property of CD4^{+}CD25^{+} T cells in allogeneic mixed lymphocyte reaction (32). One of the models adequately fitted the data, allowing separate estimates of the frequencies of proliferating and regulatory cells. In addition to the frequencies, the models provided information regarding the magnitude of the suppressor effect, allowing the detection of either one or two subtypes of regulatory cells featured by either weak or potent suppressive activity.
Materials and Methods
Singlehit Poisson model (SHPM)
The most simple model that can be fitted to limiting dilution assays (LDA) assumes that only one limiting cell of only one cell subset is necessary and sufficient for generating a positive response, defining the SHPM. According to the zero term of the Poisson distribution, the expected fraction F_{i} of negative wells is given by where f is the estimate of the cell frequency, x is the number of cells per well, i is the ith group of replicate wells.
Suppressor models describing the interaction between effector and regulatory T cells
Construction of the models is detailed in the Appendix. Briefly, we have previously proposed a model describing the interaction between regulatory T cells and effector proliferating T cells (31). This model, termed here model 1 , postulates that the inhibition of one effector proliferating T cell requires one or more regulatory cells, thus featuring a subpopulation of regulatory cells with intermediate or weak suppressive activity. These regulatory cells are termed T_{reg1} cells, in our terminology used throughout this text. The fraction F_{i} of negative wells is given by where n is the number of effector, proliferating T cells per culture (termed T_{E} cells), n ≥ 1, f_{e} is the frequency of T_{E} cells, f_{r1} is the frequency of T_{reg1} cells, Φ is a stoichiometric parameter that is defined as the minimum ratio at which T_{reg1} and T_{E} cells must be simultaneously present in the same well to ensure that this well will be negative for growth: defined for Φ ≥ 1. It must be emphasized that this parameter is not intended to provide accurate information regarding the suppressor effect of the regulatory cells on a per cell basis, because its magnitude may heavily depend on several technical parameters, like the cell density in cultures, the type or the surface of plate cultures (roundbottom wells vs flatbottom wells, 96 384well plate assay), which ultimately may influence the probability of cell contact between regulatory and effector cells. Rather, the stoichiometric parameter allows for broad discrimination between weak (large values of Φ) or intermediate (small values of Φ) suppressive activity of T_{reg1} cells. The second model, termed model 2, posits that one regulatory cell is able to inhibit one or more T_{E} cells, thus featuring a subpopulation of regulatory cells (termed T_{reg2} cells) with intermediate or potent suppressive activity. The fraction F_{i} of negative wells is written as where f_{r2} is the frequency of T_{reg2} cells, k ≥ 0. Ψ is a stoichiometric parameter that is defined as the maximum ratio at which T_{reg2} and T_{E} cells must be simultaneously present in the same well to ensure that this well will be negative for growth: defined for Ψ ≥ 1. Again, the value of Ψ allows for broad discrimination between intermediate (small values of Ψ) or potent (large values of Ψ) suppressive activity of T_{reg2} cells. The third model, termed model 3, considers that both T_{reg1} and T_{reg2} cells are present within the unfractionated cell population under study. Therefore, the fraction F_{i} of negative wells can be written as
with n ≥ 1, k ≥ 0, j ≥ 1, r ≥ 2.
Modelfitting procedure
Suppressor models 1, 2, 3 were fitted to experimental data by using the quasiNewton method to maximize the likelihood of the data (33). Briefly, the probability p(r_{i}) that r_{i} out of N_{i} wells are negative is given by where F_{i} is the theoretical fraction of negative wells given by Equations 2, 4, and 6, and r_{i} is the experimental (observed) number of negative wells. The likelihood function L for G binomial grouped observations (i.e., the number of cell dilutions performed in each LD experiment) is given by The maximum likelihood estimators of the frequency parameters f_{e}, f_{r1}, and f_{r2} are the values fê , fr1̂ , and fr2̂ which maximize where ln(L) is the natural logarithm of the likelihood function L. Notice that we used a constrained optimization procedure, imposing fixed values for the stoichiometric parameters Φ and Ψ, which can take on only integer values.
Determination of the standard error of each of the frequency parameters fê , fr1̂ , and fr2̂ with Wald statistic (34), and goodnessoffit tests evaluating the adequacy of models 1, 2, and 3 to the experimental data (33) are detailed in the Appendix.
Results
The LD data set
As an illustrative example, the modeling approach was performed on experimental LD data previously published in the literature (32). In this study, CD4^{+}CD25^{+} and CD4^{+}CD25^{−} T cells were purified from peripheral blood of healthy volunteers. Then, the cells of each subpopulation were serially diluted and cocultured with allogeneic irradiated CD3depleted stimulator cells in 96well plates for 5 days. Six (CD4^{+}CD25^{+} cells) or seven (CD4^{+}CD25^{−} cells) cell dilutions were performed, with 24 replicate wells at each cell dilution. To discriminate between positive/negative wells, cell proliferation was estimated by [^{3}H]thymidine incorporation at the end of the culture. The explicit values corresponding to the results depicted in Fig. 6 of Ng et al. (32) are given in Table I⇓. The resulting graphs according to a semilog plot representation are depicted in Fig. 1⇓. We first examined the CD4^{+}CD25^{−} LD titration curve. Although the fraction of negative cultures took the same value twice (at 156 and 313 cells per well), it appeared to be a linear relationship between the logarithm of the fraction of negative wells and the cell input. This general shape favored the singlehit kinetics hypothesis. This was checked by a generalized linear modelbased statistical slope test that we have tailored to the statistical analysis of LDA (35). Briefly, the equation of the generalized linear loglog regression line, in the form of Y_{i} = α + βX_{i} was Y_{i} = −5.874 + 0.909X_{i}. The equations of the two boundary regression lines, corresponding to the lower and upper limits of the 95% confidence interval for the slope β, were Y_{i} = −5.874 + 0.578X_{i}, and Y_{i} = −5.874 + 1.241X_{i}, evidencing that the slope β of the loglog regression line fitted to the data was compatible with 1, and ensuring that the experimental data conformed to the SHPM hypothesis. The frequency of proliferating T_{E} cells was calculated to be fê = 1/623 as estimated by the maximum likelihood method according to the SHPM hypothesis (36). We next considered the CD4^{+}CD25^{+} LD titration curve. Overall, the titration curve clearly deviates from singlehit kinetics and can be divided into three distinct phases. At low cell doses (313 and 625 cells per well), the fraction of negative wells progressively decreases and the curve exhibits approximately a linear shape (phase 1); at further increases in cell doses, the fraction of negative wells stagnates (1250 cells per well) then increases (2500 cells per well), producing a “hump” (phase 2). Then, the hump is overcome, and the fraction of negative wells decreases at 5000 cells per well (phase 3). Therefore, after the initial linear phase, the graph breaks into two regions of curvature: concave down (phase 2) and concave up (phase 3).
Fitting the suppressor models to the CD4^{+}CD25^{+} limiting dilution assay
The fitted regression lines using the bestfit frequency parameters fê , fr1̂ , and fr2̂ obtained for the three suppressor models are shown in Fig. 2⇓ along with experimental data points (x_{i}, r_{i}/N_{i}) superimposed on the graph. It can be seen that the regression lines differ primarily in the way they deal with the shape of the experimental curve. It appears that model 1fitted regression line did not fit the data well, failing to account for the three distinct phases featuring the experimental curve presented in Fig. 1⇑. By contrast, both fitted curves generated from models 2 and 3 were seemingly equivalent, yielding a good description of the experimental curve, with an initial linear phase, followed by concave down and concave up phases. Table II⇓ shows the maximum likelihood frequency estimates fê , fr1̂ , and fr2̂ , their standard errors, and Wald statistic Z scores from each of the three models. The deviance (D) and the ordinary Pearson’s χ^{2} statistic were used as goodnessoffit statistical tests to evaluate the adequacy of each model to the data. Model 1, which assumes the presence of a subpopulation of regulatory cells with intermediate or weak suppressive activity (T_{reg1} cells), suggested a high fr1̂ of 1/42, a value 14.08fold greater than the frequency of proliferating T_{E} cells, fê = 1/585. The best value of the stoichiometric parameter was Φ = 20, indicating that these regulatory cells displayed weak suppressive activity. Both fê and fr1̂ appeared to significantly contribute to the model, as evidenced by their respective Wald statistic Z scores. However, as expected from the graphical representation of the fitted regression line, model 1 poorly fits the data, with p(χ^{2}_{D} > 10.66) = 0.0137 and p(χ^{2} > 9.957) = 0.0189. Model 2, based on the presence of a subpopulation of regulatory cells with intermediate or potent suppressive activity (T_{reg2} cells), gave a low fr2̂ of 1/6258, a value 12.24fold lower than the frequency of proliferating T_{E} cells, fê = 1/511. The best value of the stoichiometric parameter was Ψ = 7, indicating that these regulatory cells displayed potent suppressive activity. Both fê and fr2̂ appeared to significantly contribute to the model, as indicated by their respective Z scores. Goodnessoffit tests gave p(χ^{2}_{D} > 7.108) = 0.0685 and p(χ^{2} > 6.694) = 0.0823, and therefore the hypothesis that this model acceptably fits the data is tenable at the 5% significance level. Model 3, based on the simultaneous presence of both T_{reg1} and T_{reg2} cell subsets, gave a very low fr1̂ of 1/10033 concomitantly with a low fr2̂ of 1/6258. It must be noticed that fr1̂ appeared to be noncontributory to the model, as indicated by the Z score, with p(z_{fr1̂} > 0.00565) = 0.995. Best values of the stoichiometric parameters were Φ = 6 and Ψ = 7. The frequency of proliferating T_{E} cells was fê = 1/511. The fact that fê and fr2̂ were identical to fê and fr2̂ of model 2, and the observation of a very low fr1̂ that did not significantly contribute to the fraction of negative wells, form the reasons why the curves, as well as the results of maximized log likelihoods, deviance, and χ^{2} from goodnessoffit tests, are identical for both models 2 and 3. However, model 3 includes only two degrees of freedom (g = 3) instead of three (g = 2) for model 2, because of the inclusion of fr1̂ as an additional parameter. It was found p(χ^{2}_{D} > 7.108) = 0.0286 and p(χ^{2} > 6.694) = 0.0351 for model 3, and therefore the hypothesis that this model acceptably fits the data was rejected. In conclusion, among these three models, we retain model 2 as the unique model that acceptably fits the CD4^{+}CD25^{+} LD data; 95% CI for fê was 1/837 − 1/368 and 95% CI for fr2̂ was 1/15889 − 1/3894.
It was of interest to determine whether extension of model 2 with an additional parameter could significantly improve the fit to the CD4^{+}CD25^{+} LD data. To this end, we took into account the possibility that a fraction of the titrated cells might proliferate, escaping from the suppressive activity of the regulatory cells. These cells refractory to growth inhibition might include a minor fraction of proliferating T_{E} cells, but also a part of the regulatory cells themselves. The possibility that strongly activated, proliferating T_{E} cells can be refractory to the inhibitory effect of CD4^{+}CD25^{+} regulatory cells has been reported (18). The presence of these suppressorresistant cells can be modeled by simply multiplying Equation 4 by the term exp(−f_{res}x_{i}), which is the zero term of the Poisson distribution with parameter f_{res}, the sum of the frequencies of suppressorresistant cells that might originate from proliferating T_{E} and regulatory cell subtypes. Therefore, the fraction F_{i} of negative wells can be written as where f_{res} is the frequency of proliferating, suppressorresistant cells, with f_{er} and f_{sr} denoting the frequencies of proliferating, suppressorresistant cells among proliferating T_{E} and regulatory cells, respectively. Maximum likelihood estimate of f_{res} under this new model nested within model 2 gave freŝ = 0. Therefore, the nested model reduced to model 2, invalidating the hypothesis of suppressorresistant cells among the titrated CD4^{+}CD25^{+} cells.
Discussion
Several subsets of regulatory cells with distinct phenotypes and distinct mechanisms of action have now been identified among CD4^{+} T cells. These include type 1 (Tr1) regulatory cells, which secrete high levels of IL10 and low to moderate levels of TGFβ; Th3/Tr2 regulatory cells, which primarily secrete TGFβ; and CD25^{+} regulatory T cells, which inhibit immune responses primarily through cellcell contact (1, 2, 3, 4, 5, 6, 7, 8, 9). Although their suppressor effect, on a per cell basis, is not accurately known, it is usually considered that the CD4^{+}CD25^{+} regulatory cells display a remarkably potent suppressive effect on proliferation and cytokine production by both CD4^{+} and CD8^{+} T cells (37). In vitro, the most widespread protocol consists of coculturing CD4^{+}CD25^{−} responder T cells and CD4^{+}CD25^{+} T cells, in the presence of soluble or platebound antiCD3 Abs and irradiated APC. In multiple experiments of this type, complete suppression is usually observed at a ratio of suppressors/responders of 1:1, and significant suppression (>70%) can persist at a ratio of 1:4 (37, 38, 39, 40) and even at 1:10, as demonstrated recently with tumorspecific regulatory CD4^{+}CD25^{+} T cell clones (41). It is generally considered that Tr1 cells have a less potent suppressive activity than CD25^{+} regulatory cells, mainly because they mediate their action by soluble suppressive cytokines (IL4, IL10, TGFβ) instead of cell contact. Thus, in our models, T_{reg1} cells may be equivalent to Tr1 cells, and T_{reg2} cells to CD4^{+}CD25^{+} regulatory cells. The current modeling study of the CD4^{+}CD25^{+} LD data has favored the hypothesis of a single regulatory cell subset with potent suppressive activity (Ψ = 7), and therefore these results are in full agreement with the literature. More surprising, it was observed that the frequency (f_{e}) of proliferating cells among the CD4^{+}CD25^{+} T cells was comparable to that observed among the CD4^{+}CD25^{−} T cells (1/511 versus 1/623, respectively). This is an unexpected result, since it is well established that activated naturally occurring CD4^{+}CD25^{+} T cells cannot proliferate in an autonomous fashion, as a direct consequence of their inefficiency to produce their own autocrine growth factors like IL2 (42, 43). It must be noticed that the CD4^{+}CD25^{+} T cells were not extensively purified in the current LD experiment, with a purity ranging from 80 to 90% (32). Thus, one plausible explanation is that the observed frequency of proliferating cells may be in part attributable to the contaminating alloantigenactivated CD4^{+}CD25^{−} T cells. On the other hand, the hypothesis that a fraction of naturally occurring CD4^{+}CD25^{+} T cells may be actually effector proliferating cells rather than regulatory cells has yet to be evoked. In this view, a study performed at the clonal level (44) has highlighted the functional heterogeneity of this cell population. Although derived from extensively purified CD4^{+}CD25^{high} polyclonal T cells, which were claimed to represent a pure subpopulation of regulatory cells (18, 45, 46), a notable fraction of CD4^{+}CD25^{+} T cell clones (at least 20%) were not anergic, produced IL2, and were unable to suppress the proliferative response of autologous CD4^{+}CD25^{−} T cells. LD experiments performed with highly purified CD4^{+}CD25^{+} T cells should be of great help to resolve the question of the functional heterogeneity of this cell subset. Finally, the most interesting point of our current study relates to the remarkably low frequency of T_{reg2} cells (1/6258), which were functionally engaged in the allogeneic response. This result suggests that regulatory cells actually implicated in a given immune response are scarce, representing an extremely low percentage of the CD25^{+} cell subset, and therefore emphasizing that the regulatory cell compartment cannot be reliably estimated with this marker. In addition, these findings suggest that, at least in this LD experiment, the activation of regulatory cells was Agspecific, and that they did not convey suppressor activity to CD4^{+}CD25^{+} lymphocytes other than those directed against the initial alloantigens. In the literature, the conversion of CD4^{+}CD25^{−} effector, proliferating cells to regulatory cells under the influence of naturally occurring CD4^{+}CD25^{+} regulatory cells has been reported, and this process is known as “infectious tolerance” (47, 48) or “bystander immune suppression” (49). This mechanism has not been demonstrated within the natural CD4^{+}CD25^{+} regulatory cells themselves.
Numerous studies in the literature have reported examples of LD experiments that do not conform to singlehit kinetics. These nonsinglehit, multihit kinetics curves can arise in a variety of experimental systems including allogeneic or autologous mixed lymphocyte reactions, Ag or mitogeninduced T cell proliferation, T cellmediated cytotoxicity assays (29, 32, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60). The resulting LD titration curves exhibit variable shapes. Typically, at the lowest cell doses, the logarithm of the fraction of negative wells progressively decreases, appearing to be linearly related to the cell dose. Then, as the cell number increases, this initial linear segment with a negative slope is followed by a kink, and the fraction of negative wells increases, generating a segment with a positive slope. In a number of examples, the fraction of negative wells can increase up to 100%. Such LD titration curves have been termed “Vshaped curves” by Lefkovits and Waldmann (29). Another possibility is that the fraction of negative wells stagnates after the kink, generating a plateau (leveling off curves). More complicated curve shapes have been frequently described, with three distinct segments (29, 32, 52, 53, 54, 57, 58, 59, 60, 61): the initial Vshaped segment is followed by a second kink, then the fraction of negative wells decreases as the cell number is increased, generating an ultimate segment with a negative slope.
More than 20 years ago, the nonlinear appearance of LD titration curves was interpreted as resulting from interactions among effector and regulatory cell subtypes (50, 51, 52, 53). However, very few studies in the literature have been committed to model such limiting dilution data, in search of regulatory cells mixed with effector cells. Lefkovits and colleagues (29, 60) have proposed a number of mathematical models yielding nonlinear titration curves but no attempt was made to fit these models to experimental data with suitable mathematical methods. Fey et al. (52, 61) proposed a model based on the presence of one subpopulation of regulatory cells with weak suppressive activity (thus equivalent to T_{reg1} cells with frequency f_{r1}, according to our notations) and described by This model resembles model 1, Equation 2, because the term is equivalent to the second term on the righthand side of Equation 2 except that a function A(n) has been introduced. The function A(n) was defined as with y(n) = or y(n) = Γn, or y(n) = Γ (52, 53, 54). This function means that the minimum number A(n) of T_{reg1} cells that must be simultaneously present with n proliferating cells in the same well to suppress their growth, is not directly proportional to n as described in model 1, which assumes T_{reg1} = ΦT_{e} (corresponding to Equation 3). In a variety of LD experiments, Γ was found to be in order from 5 to 20, with the value of a ranging from 5 to 10, assuming that there was rather a large number of regulatory cells required for suppression of one proliferating cell. Moreover, the results of the modeling studies suggested a notable excess of T_{reg1} cells relatively to proliferating T_{E} cells, with (f_{r1} /f_{e}) > 20 (52, 53, 54). However, these data must be currently considered cautiously. Equation 8 was not solved numerically by optimization methods, and “best” values of all parameter estimates, f_{e}, f_{r11}, a, and Γ, were empirically determined. Moreover, it was not presented any statistical evaluation of the goodnessoffit of the model to the LD data.
More recently, another model has been proposed by Dozmorov et al. (57). This model postulated that a fixed number of regulatory cells, termed a, was able to inhibit the growth of an unlimited number of proliferating T_{E} cells, and that these regulatory cells proliferated in turn when their number reached a fixed value in wells. Thus, this model is based on the presence of one subpopulation of regulatory cells with potent suppressive activity (equivalent to T_{reg2} cells with frequency f_{r2}) and is written as according to our notations; b is the threshold number per well of regulatory cells at which these cell can proliferate, with b > a. The twomodel parameters a and b were properly estimated from the data by a numerical optimization procedure. This model leads to “zigzag” LD titration curves, getting their name due to the shape of the curve, exhibiting three distinct phases. Typically, at the lowest cell doses, the fraction of negative wells progressively decreases according to singlehit kinetics, then there is a first kink and the fraction of negative wells increases at intermediate cell doses. Finally, after a second kink, the fraction of negative wells ultimately decreases at the highest cell doses. Although this model can yield a good fit to LD data (57, 58), one major difficulty arises in attempt to explain the mechanism(s) by which a given fixed number of regulatory cells might be able to suppress the growth of an unlimited number of proliferating T_{E} cells, and therefore the biological relevance of this too simplistic model is questionable. Moreover, the model nested within model 2 presented in our current study (Equation 7) especially fits well such zigzag curves, as illustrated with a set of theoretical curves (data not shown). Over the Dozmorov model, this model offers the advantage to make a simple, and reasonable assumption of a subset of proliferating suppressorresistant cells.
Although model 3 did not appear to fit well to the CD4^{+}CD25^{+} LD data, it must be noticed that, apart from allowing the estimation of the frequencies of two distinct subpopulations of regulatory cells, it theoretically offers the opportunity to evaluate their possible synergistic effect on the inhibition of T cell growth. Let us consider the fourth term on the righthand side of Equation 6 This term describes the possibility that a well will be negative for growth if it contains both T_{reg1} and T_{reg2} cells, in the situation where neither T_{reg1} nor T_{reg2} cells reach the critical number of cells required to inhibit all T_{E} cells per se (see Appendix). A well will be negative for growth if it contains a sufficient number, m, of T_{reg1} cells, able to suppress the growth of the remaining number of T_{E} cells, which are not inhibited by T_{reg2} cells. The condition relative to m is where Φ is the minimum number of T_{reg1} cells required to inhibit one T_{E} cell, Equation 3. It can be easily seen that a synergistic mechanism between T_{reg1} and T_{reg2} cells should be approached by simply decreasing the value of Φ in this term, thus decreasing the number m of T_{reg1} cells required to inhibit the growth of the remaining number of T_{E} cells, which are not inhibited by T_{reg2} cells. The condition relative to m becomes with Φ_{s} < Φ. Several values of Φ_{s} should be tested, in search of the best fit of model 3 to the data. Although the present study was illustrated in the context of CD4^{+}CD25^{+} regulatory cells, the models developed here are primarily intended to explore the regulatory cell compartment in unfractionated cell population, irrespective of their lineage (Tr1 cells, Th3 cells, naturally occurring Foxp3expressing CD4^{+}CD25^{+}, or TGFβconverted Foxp3expressing CD4^{+}CD25^{+} regulatory cells derived from naive Foxp3negative CD4^{+}CD25^{−} T cells). Apart from the accurate quantitation of regulatory cells, this approach allows for detection of either one or two separate regulatory cell subsets, based on the magnitude of their suppressor activity. In conclusion, although mathematical modeling has yielded some quantitative results that have improved our understanding of the immune system (62, 63), this approach is considered as still in its infancy (64, 65). The recent successes of models of HIV1 and T lymphocyte dynamics in HIV1 seropositive patients (65) provide a strong motivation for pursuing mathematical approaches, applied to model the functional interactions between effector and regulatory cells in immune responses. Limiting dilution data modeling should be a tool of choice to study these interactions. Although more elaborate mathematical models that parallel the rapid increasing knowledge regarding regulatory cells will be conceivable in the future, our aim was to define generic equations providing a basis for modeling limiting dilution data.
Disclosure
The authors declare no financial interest.
Appendix 1
Construction of the suppressor models describing the interaction between effector and regulatory T cells and applied to a limiting dilution assay according to a Poisson process
Model 1
The first assumption that a well will be scored negative for growth is that it did not receive any T_{E} cell. The probability of this event is given by Next, as the second assumption, a well will be negative for growth if it received simultaneously, n T_{E} cells and at least (n × Φ) T_{reg1} cells, i.e., exactly (n × Φ), or more than (n × Φ), T_{reg1} cells. Thus, a first possibility that a well will be scored negative for growth is that it received one T_{E} cell and at least Φ T_{reg1} cells (i.e., exactly Φ, or >Φ, T_{reg1} cells). The probability of this event is given by A second possibility is that it received two T_{E} cells and at least (2 × Φ) T_{reg1} cells, i.e., exactly (2 × Φ), or >(2 × Φ), T_{reg1} cells. The probability of this event is given by and so on. Thus, the total fraction of negative wells is given by the sum of probabilities over all n to obtain a negative well, and can be adequately described by the general formula given by Equation 2.
Model 2
Apart from the assumption that a well did not receive any T_{E} cell given by the term exp(−f_{e}x_{i}), the assumption that a well will be negative for growth is that it received up to (Ψ × k) T_{E} cells and simultaneously, at least k T_{reg2} cells (i.e., exactly k or >k T_{reg2} cells). Thus, in addition to exp(−f_{e}x_{i}), a first possibility to obtain a negative well is that it received 1, 2, …, Ψ T_{E} cells (i.e., from one to Ψ T_{E} cells) and at least, one T_{reg2} cell (i.e., exactly 1 or >1 T_{reg2} cells). The probability of this event is given by A second possibility is that it received (Ψ + 1), (Ψ + 2), …, (Ψ + Ψ) T_{E} cells, i.e., from (Ψ + 1) to (2 × Ψ) T_{E} cells and at least, two T_{reg2} cells (i.e., exactly 2 or >2 T_{reg2} cells). The probability of this event is given by A third possibility is that it received (2 × Ψ + 1), (2 × Ψ + 2), …, (2 × Ψ + Ψ) T_{E} cells, i.e., from (2 × Ψ + 1) to (3 × Ψ) T_{E} cells and, at least, three T_{reg2} cells (i.e., exactly 3 or >3, T_{reg2} cells). The probability of this event is given by and so on. Thus, the total fraction of negative wells is given by the sum of probabilities over all k to obtain a negative well and can be adequately described by the general formula given in Equation 4.
Model 3
In addition to the term exp(−f_{e}x_{i}), this model posits that a well will be negative for growth according to three assumptions: the presence of a sufficient number of T_{reg1} cells as described by model 1 and irrespective of the number of T_{reg2} cells (assumption A); the presence of a sufficient number of T_{reg2} cells as described by model 2 and irrespective of the number of T_{reg1} cells (assumption B). The third hypothesis considered by model 3 is that a well contains both T_{reg1} and T_{reg2} cells but that neither T_{reg1} nor T_{reg2} cells reach the critical number of cells required to inhibit all T_{E} cells per se but that together they can; therefore, the well will be negative for growth (assumption C).
We first consider assumptions A and B. The probability of the assumption A, called event A, is given by the second term on the righthand side of Equation 2, i.e., the suppressor term of model 1 The probability of the assumption B, called event B, is given by the second term on the righthand side of Equation 4, i.e., the suppressor term of model 2 A and B are independent events; that is, A does not change calculation of the probability of the event B and reciprocally. However, it must be noted that A and B are not mutually exclusive, that is, have outcome in common: a well will be negative if it contains n T_{E} cells and a sufficient number of T_{reg1} cells, but the same well may simultaneously contain a sufficient number of T_{reg2} cells to inhibit per se the growth of the n T_{E} cells. The intersection of these two events, A ∩ B, is the event that both A and B occur in the same well. Thus, the correct probability that a well will be negative for growth is given by and according to multiplication rule for independent events, it can be written and therefore model 3 can now be written as
Determination of p(A) × p(B)
A first possibility corresponding to p(A) ∩ p(B) is that a well had received, 1, 2, …, Ψ T_{E} cells, and at least one T_{reg2} cell (i.e., p(B)) and at least (n × Φ) T_{reg1} cells (i.e., p(A)), where n is the number of T_{E} cells, n = 1, 2, …, Ψ. The probability of this event is given by A second possibility corresponding to p(A) ∩ p(B) is that a well had received (Ψ + 1), (Ψ + 2), …, (Ψ + Ψ) T_{E} cells, and at least two T_{reg2} cells (i.e., p(B)), and at least (n × Φ) T_{reg1} cells (i.e., p(A)), where n is the number of T_{E} cells, n = (Ψ + 1), (Ψ + 2), …, (Ψ + Ψ). The probability of this event is given by A third possibility corresponding to p(A) ∩ p(B) is that a well had received (2 × Ψ + 1), (2 × Ψ + 2), …, (2 × Ψ + Ψ) T_{E} cells, and at least three T_{reg2} cells (i.e., p(B)) and at least (n × Φ). T_{reg1} cells (i.e., p(A)), where n is the number of T_{E} cells, n = (2 × Ψ + 1), (2 × Ψ + 2), …, (2 × Ψ + Ψ). The probability of this event is given by and so on. Thus the general term corresponding to p(A) ∩ p(B) is the sum of all these probabilities over all k and is given by We now examine assumption C, i.e., p(C), to obtain a negative well. To construct the general term that adequately describes this assumption, let us consider the possibility that a well contains n = (Ψ + 1), (Ψ + 2), (Ψ + 3), …, (Ψ + r) T_{E} cells, and exactly one T_{reg2} cell. The probability of this event can be written as where r ≥ 2, r taking a value in the set Ω = {2, 3, …, R}. Next, the crucial point to consider is that the number of remaining T_{E} cells that are not inhibited by the T_{reg2} cell is given by Therefore, the well will be negative for growth if it contains a sufficient number, m, of T_{reg1} cells able to suppress the growth of the remaining number of T_{E} cells that are not inhibited by the T_{reg2} cell. This number, m, is defined as where Φ is the minimum number of T_{reg1} cells required to inhibit one T_{E} cell, Equation 3. Thus, the probability to have a negative well is given by Now, according to the hypothesis of model 3, recall that, like T_{reg2} cells, the number of T_{reg1} cells must not be sufficient to inhibit per se the growth of all T_{E} cells. Thus, there is a constraint on m, and finally the probability to have a negative well is given by A second possibility is that a well contains n = (2 × Ψ + 1), (2 × Ψ + 2), (2 × Ψ + 3), …, (2 × Ψ + r) T_{E} cells and exactly two T_{reg2} cells. The probability of this event can be written as with r ≥ 2, and the number of remaining T_{E} cells that are not inhibited by the two T_{reg2} cells is given by (n − 2Ψ). Therefore, the well will be negative for growth if it contains a sufficient number, m, of T_{reg1} cells able to suppress the growth of the remaining number of T_{E} cells that are not inhibited by the two T_{reg2} cells. The condition on m is and finally the probability to have a negative well is given by and so on. Thus the general term corresponding to the probability to have a negative well according to the assumption C of model 3 is given by the sum of probabilities over all j and is given by where j is the number of T_{reg2} cells, and finally the general formula describing model 3 is given by Equation 6. In practice, the starting value of r is 2, then incremented until an appropriate value is set, that is, the ultimate value of r that causes a significant contribution to the probability considered above.
Determination of the standard error of each of the frequency parameters fê , fr1̂ , and fr2̂ and Wald statistic
Calculation of the standard error of each of the frequency parameters fê , fr1̂ , and fr2̂ needs to construct the matrix of the second derivatives of the likelihood function L with respect to its parameters (34 ). This matrix is commonly referred to as the information matrix or Hessian matrix and contains the negative second partial derivatives of ln(L) with respect to each of the frequency parameters fê , fr1̂ , and fr2̂ . For models 1 and 2, which are twoparameter ( fê and fr1̂ ) models, the Hessian matrix is written as and for model (3), which is a threeparameter ( fê , f̂_{r1}, f̂_{r2}) model, the Hessian matrix is written as The matrix inverse of H, named Var̂ (f̂), provides the estimated variancecovariance matrix of the maximum likelihood estimates fê , fr1̂ , and fr2̂ . The estimated SE of each of the frequency estimates fê , fr1̂ , and fr2̂ is the square root of each diagonal element of (H)^{−1}, and the offdiagonal elements of (H)^{−1}, which are not reported below, are the covariances between pairs of frequency estimates. For the twoparameter models 1 and 2 and, for the threeparameter model 3 The Wald test is a significance test of whether the frequency estimate f̂ of interest is significantly different from 0. With nonnull SE of f̂, the test statistic is and has an approximate standard normal distribution when f̂ = 0. When f̂ ≠ 0, an approximate 95% confidence interval for f̂ can be calculated as
Goodnessoffit tests evaluating the adequacy of models 1, 2, 3 to the experimental data (33 )
Let r̂_{i} be the fitted number of negative wells under the current model, r̂_{i} = N_{i}F̂_{i}, where F̂_{i} is the fitted fraction of negative wells. The deviance can be written as comparing the observations r_{i} with their corresponding fitted values rî under the current model. The Pearson χ^{2} statistic is Both the deviance and the Pearson χ^{2} statistic have the same asymptotic χ^{2} distribution, with (G − g) df, where g is the number of parameters estimated from the data by the modelfitting procedure. All calculations were performed with Mathematica 4.2 software (Wolfram Research) and Global Optimization 4.2 software written for Mathematica (Loehle Enterprises).
Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

↵1 This work was supported by Grant 3218 from Association pour la Recherche sur le Cancer.

2 Address correspondence and repinrt requests to Dr. Thierry Bonnefoix, Institut National de la Santé de la Recherche Médicale U353, Institut Albert Bonniot, RondPoint de las Chantourne, 38706 Grenoble, France. Email address: Thierry.Bonnefoix{at}ujfgrenoble.fr

↵3 Abbreviations used in this paper: GITR, glucocorticoid TNF receptor familyrelated gene; LD, limiting dilution; LDA, LD assay(s); SHPM, singlehit Poisson model; T_{E}, effector T cell(s); T_{reg1} and T_{reg2}, types 1 and 2 regulatory cells.
 Received May 18, 2004.
 Accepted January 6, 2005.
 Copyright © 2005 by The American Association of Immunologists