Supplementary Materials
Here, details are given regarding the mathematical models and methodologies that have been used to arrive at the conclusions described in the main text.
The data
All data used to perform calculations in this study were taken from a previous publication (Roeder et al. 2006). These data document the time course of BCR-ABL1/ABL1 levels during imatinib treatment (400mg/day) in 69 CML patients that were newly diagnosed, recruited from the German cohort of the IRIS study. The patients were followed for 5 years. Because the data for the individual patients are provided in (Roeder et al. 2006), they are not re-displayed here.
Quantifying growth and decline slopes in the individual patient data
Overall, three types of slopes can be found in the individual patient data: the slope of a first and fast phase of decline, the slope a second and slower phase of decline, and a growth slope. In most patients, all three slopes are found. Some patients can have a subset of these. Before the first and fast phase of decline, many patients show a “shoulder phase” during which no significant decline in BCR-ABL levels is observed, or during which the BCR-ABL levels even increase. This phase has been ignored in the following procedures. To quantify the three slopes mentioned above, a three phase exponential growth/decline model was used, given by F = a1 exp(b1t) + a2 exp(b2t) + a3 exp(b3t), where bi represent the slopes of exponential growth/decline. Using non-linear least squares regression and standard software, this model was fitted to the individual patient data and the three slopes were recorded. If the CML decline dynamics were characterized by less than these three phases (e.g. only one decline phase was present), this was reflected in the results from the data fitting, and only the appropriate slopes were recorded. In order to correlate the different slopes with each other, only the patients that were characterized by the presence of the appropriate decline phases were used.
Mathematical model of immune – CML dynamics
A mathematical model was used that was previously published in the context of immune-impairing virus infections, such as HIV (Komarova et al. 2003). The model contained two variables: a growing virus population that has the ability to impair the immune system, and a specific immune response. In the present context, we are interested in CML and thus consider a population of growing CML cells that have the potential to impair immunity, and an immune response that reacts against CML antigens. The equations are the same and are given by:
dx/dt = rx(1-x/k) – dx – pxz;
dz/dt = cxz/[ (z+h)(x+e) ] – qxz – bz.
The CML population is denoted by x, while the immune cell population is denoted by z. The CML population grows logistically. That is, at low numbers of cells, growth is exponential while growth slows down as the number of cells increases. The tumor cells die with a rate d and become removed by the immune response with a rate p. The immune response expands upon antigenic stimulation by tumor cells with a rate c. Immune expansion saturates if the number of tumor cells and the number of immune cells are high. The tumor cells impair immunity with a rate q. Finally, immune cells decay in the absence of antigenic stimulation with a rate b. Data suggest that the relevant immune response observed during CML therapy is a T cell response consisting of both CD4+ and CD8+ T cells (Chen et al. 2008; Kim et al. 2008; Wang et al. 2005). T cell responses can be modeled in a variety of ways (Antia et al. 2005; De Boer & Perelson 1998; Wodarz & Thomsen 2005). However, the model presented above is merely an example of a class of models that share common mathematical properties (Komarova et al. 2003). Therefore, the modeling approach taken here does not depend on the particular mathematical formulations written down in the model above. Results derived from this model remain robust as long as the following general assumptions remain true (Komarova et al. 2003): The tumor cells can both stimulate and impair their specific immune responses; Immune expansion is negative if the number of cancer cells is very low due to lack of significant antigenic stimulation; Immune expansion is also negative if the number of tumor cells is relatively high because of overwhelming immune impairment; Immune expansion is positive for an intermediate number of cancer cells because the degree of immune stimulation outweighs the degree of immune impairment. Mathematical properties of this class of models and aspects of robustness have been summarized briefly in the text and are described in detail in (Komarova et al. 2003).
For the current study, the above described model was adapted to include CML therapy and the presence of both drug-sensitive and drug-resistant mutants. Denoting drug sensitive tumor cells by x and drug resistant tumor cells by y, the model is given by the following ordinary differential equations.
dx/dt = -dx – pxz
dy/dt = ry – pyz
dz/dt = c(x+y)z/[ (z+h)(x+y+e) ] – q(x+y)z – bz.
Because drug-sensitive cells, x, are susceptible to therapy, this population of cells is assumed to decline during treatment with a rate d. As before, immune responses contribute to cell death with a rate p. On the other hand, drug-resistant cells, y, are not susceptible to therapy and are thus assumed to expand exponentially during therapy with a rate r. Because the initial growth phase of the resistant CML cells is interesting in the current context, we ignore growth saturation at high numbers of cells for simplicity. As with drug-sensitive cells, drug-resistant cells are removed by immune responses with a rate p. The equation for immune responses is the same as in the simple model described previously. However, both antigenic stimulation and immune impairment are now proportional to the number of drug-sensitive and drug-resistant cells, x+y.
Non-linear least square fitting of model to data
The mathematical model was fitted to the individual data sets using non-linear least square regression. This was performed with a software package, Berkeley Madonna (http://www.berkeleymadonna.com/). The model can fit a diverse set of clinical profiles and is thus consistent with observation. The biological meaning of the parameter values that were obtained from the fit is, however, unclear. The fitting procedure finds a minimum at which the discrepancy between the data and the predicted curve is lowest. Usually, however, there is not just one such minimum, but a collection of several local minima. In other words, several parameter combinations can exist that result in a comparably good fit of the model to the data. Therefore, the parameter values that are obtained from the fit should not be interpreted as biologically “true” parameter measures without further experimental work. This is a limitation inherent in the method. However, the fitting does show that the model can describe the observed data sets, and this should encourage further work that tries to directly measure individual parameters with experimental methods.
Model predictions about correlations between different slopes
According to the model, the slope of the fast and slow phase of decline are influenced by the temporary rise in immunity as therapy is started. In addition, the slope of slow phase of decline and the slope of the tumor re-growth are influenced by the presence of drug-resistant mutants in the tumor cell population. To examine the correlations predicted by the model, computer simulations were run many times, randomly varying the rate of immune proliferation (parameter c) and the replication rate of the resistant mutants (parameter r).
The slopes were quantified as follows in the computer simulation. The slopes are not constant during the various phases of decline and growth, but vary continuously. This makes the terms “fast phase” and “slow phase” somewhat arbitrary. Nevertheless, this classification is useful in the current context. For the fast phase of decline, the slope first increases and then starts to decrease. The fastest slope of decline was recorded in the computer simulation. For the slope of the second and slower phase of decline, the average slope between the fastest slope and the time point when the tumor cell population began to re-grow was calculated. The slope of the eventual rise was simply given by the growth rate of the drug-resistant tumor cell population, r. When the tumor began to grow in the simulation, immunity had decayed to insignificant levels and the tumor cell population initially grew with a rate r. Of course, later during the growth phase, it is possible for the immune response to expand again (depending on the parameter values), and counter the exponential growth of the cancer cell population, as illustrated in Figure 4 in the main text. However, since the analysis focused on the initial re-growth of the tumor cells, it is appropriate to just consider the initial exponential rate of increase.
In the computer simulations, parameter combinations were excluded that were irrelevant to predicting the correlations between the three slopes. These include the following cases: absence of drug-resistant mutants; failure of immune responses to rise and significantly contribute to the decline phases; high initial abundances of drug-resistance mutants such that they are the main determinants of the CML decline during the early stages of therapy; presence of essentially only a single phase of CML decline during treatment.
Using these calculations, it was found that the model with its underlying assumptions can predict the observed correlations between the slopes: there is a positive correlation between the slope of the fast and slow phase of decline, a negative correlation between the slope of the slow phase of decline and the rise, and a lack of correlation between the slope of the fast phase of decline and the rise.
However, if one of the variables is varied over a significantly wider range than the other and is thus the dominant factor driving the dynamics, different correlations can be observed. If there is minimal variation in the growth rate of the resistant mutant, the positive correlation between the slopes of the fast and the slow phase of decline is maintained (because it is determined by variation in the strength of immunity), but the negative correlation between the slopes of the slow phase of decline and the rise disappears, leading to a lack of correlation. If there is minimal variation in the strength of immunity, the growth rate of resistant mutants becomes the only significant driving force of the dynamics. This can have various complicating effects on the correlations. If there is no significant variation in the strength of immunity, variations in the slope of the fast decline phase can only be brought about by the presence of the resistant mutant. The faster it grows, there higher the early levels of the resistant mutant and the slower the fast phase of decline. Since the abundance of the resistant mutant at the beginning of therapy is likely to be small, this effect is also likely to be small. During later phases of treatment, the abundance of the resistant mutant will be higher, and can slow down the second phase of decline. As before, this can again lead to a positive correlation between the fast and slow decline slopes. On the other hand, the resistant mutant can also stimulate the immune response. A faster growth rate of the resistant mutant can then lead to more immune stimulation and thus to higher levels of immune responses, which in turn can speed up the second phase of decline. This would lead to a negative correlation between the fast and slow decline slopes, as well as to a positive correlation between the slopes of the slow decline phase and the rise. However, these are only mentioned for completeness and are likely to be irrelevant from a biological point of view. According to the theory developed here, the clinical data can be explained if there is pronounced variation both in the strength of immunity and in the growth rate of drug-resistant mutants, which gives rise to the correlations described in the main text.
Finally, while variation in the strength of the immune response has been expressed by varying the rate of immune proliferation, c, it could have also been expressed by varying the death rate of the immune cells, b. A high value of b correlates with a weaker response, while a lower value of b correlates with a stronger response. Variation in the parameter describing the degree of immune impairment, q, is not meaningful in this context because the dynamics described here assume that the degree of cancer-induced immune impairment lies within an intermediate range. If it is too strong, then an immune response cannot rise. If there is hardly any immune impairment, the immune response would control the cancer in the first place.
Antia, R., Ganusov, V. V. & Ahmed, R. 2005 The role of models in understanding CD8+ T-cell memory. Nat Rev Immunol 5, 101-11.
Chen, C. I., Maecker, H. T. & Lee, P. P. 2008 Development and dynamics of robust T-cell responses to CML under imatinib treatment. Blood 111, 5342-9.
De Boer, R. J. & Perelson, A. S. 1998 Target cell limited and immune control models of HIV infection: a comparison. J Theor Biol 190, 201-14.
Kim, P. S., Lee, P. P. & Levy, D. 2008 Dynamics and potential impact of the immune response to chronic myelogenous leukemia. PLoS Comput Biol 4, e1000095.
Komarova, N. L., Barnes, E., Klenerman, P. & Wodarz, D. 2003 Boosting immunity by antiviral drug therapy: a simple relationship among timing, efficacy, and success. Proc Natl Acad Sci U S A 100, 1855-60.
Roeder, I., Horn, M., Glauche, I., Hochhaus, A., Mueller, M. C. & Loeffler, M. 2006 Dynamic modeling of imatinib-treated chronic myeloid leukemia: functional insights and clinical implications. Nat Med 12, 1181-4.
Wang, H., Cheng, F., Cuenca, A., Horna, P., Zheng, Z., Bhalla, K. & Sotomayor, E. M. 2005 Imatinib mesylate (STI-571) enhances antigen-presenting cell function and overcomes tumor-induced CD4+ T-cell tolerance. Blood 105, 1135-43.
Wodarz, D. & Thomsen, A. R. 2005 Does programmed CTL proliferation optimize virus control? Trends Immunol 26, 305-10.
7