Why do low-level exposures to environmental toxins often elicit over-compensating responses that reduce risk to an organism? Conversely, if these responses improve health, why wait for an environmental challenge to trigger them? This paper presents a mathematical modeling framework that addresses both questions using the principle that evolution favors tissues that hedge their bets against uncertain environmental challenges. We consider a tissue composed of differentiated cells performing essential functions (e.g., lung tissue, bone marrow, etc.). The tissue seeks to maintain adequate supplies of these cells, but many of them may occasionally be killed relatively quickly by cytotoxic challenges. The tissue can “order replacements” (e.g., via cytokine network signaling) from a deeper compartment of proliferative stem cells, but there is a delivery lag because these cells must undergo maturation, amplification via successive divisions, and terminal differentiation before they can replace the killed functional cells. Therefore, a “rational” tissue maintains an inventory of relatively mature cells (e.g., the bone marrow reserve for blood cells) for quick release when needed. This reservoir is replenished by stimulating proliferation in the stem cell compartment. Normally, stem cells have a very low risk of unrepaired carcinogenic (or other) damage, due to extensive checking and repair. But when production is rushed to meet extreme demands, error rates increase. We use a mathematical model of cell inventory management to show that decision rules that effectively manage the inventory of mature cells to maintain tissue function across a wide range of unpredictable cytotoxic challenges imply that increases in average levels of cytotoxic challenges can increase average inventory levels and reduce the average error rate in stem cell production. Thus, hormesis and related nonlinearities can emerge as a natural result of cell-inventory risk management by tissues.
Hormesis has been proposed as a widespread phenomenon (Calabrese, 2005), cutting across many organisms, causal agents and mechanisms, and types and levels of biological responses, including cancer (Fukushima et al., 2005). Yet, the idea that organisms require exposures to low levels of environmental stressors to achieve optimal health may not seem intuitively very plausible. Why would evolution lead to such a result, especially for life-threatening responses such as cancer? How can reduction of cancer risks by exposures to known carcinogens be reconciled with evolutionary pressures to select organisms that are well fitted to their environments?
This paper explores the following possible answer: homeostasis, which is essential for survival in unpredictable environments, implies the possibility of hormesis and related dose-response nonlinearities for many adaptive systems and exposure patterns. According to this proposed explanation, biological systems that are able to cope with unpredictable stresses and maintain adequate functioning over a wide range of exposure conditions should be expected to exhibit non-linear responses to stresses. Some types of non-linearities (e.g., J-shaped and U-shaped dose-response relations) commonly referred to as hormesis arise naturally in such feedback control systems.
This paper develops this potential explanation of hormesis in the context of the dynamic responses of cell populations in tissues or organs exposed to cytotoxic agents. In such systems, effective feedback control is essential for replenishing cell populations depleted by cell-killing. While detailed physiologically-motivated models of cell kinetics in bone marrow and blood have previously been developed and compared to experimental and clinical data (e.g., Steinbach et al., 1980; Atkins, 2000), they include explicitly non-linear components and control laws, so that it is not surprising that they exhibit nonlinear responses to cytotoxic doses. A methodological contribution of the present paper is to show that important nonlinearities arise even in adaptive systems with linear underlying dynamics when cell populations are replenished by delayed production of new cells, as is physiologically realistic. Such delayed feedback systems exhibit many of the key qualitative nonlinearities observed in previous more complicated simulation models and in experimental and clinical data (Atkins, 2000; Cox, 2000).
METHODS: A SIMPLIFIED MODEL OF HOMEOSTASIS
A simple, concrete model of an adaptive system is very useful for examining how nonlinear responses can arise in systems with delayed feedback control. To this end, we suppress details and greatly simplify a previously published, physiologically-based model of hematotoxicity dynamics in humans, mice, and dogs (Cox, 2000). The role of delayed feedback control in producing nonlinear responses is made as clear as possible by modeling the normal production of cells as a linear timeinvariant dynamic system (represented by a system of linear ordinary differential equations with constant coefficients) thus suppressing nonlinearities such as s-shaped curves in rates of cell production in the original model. Figure 1 shows the structure of the simplified model.
Structure of a simplified model of cell production and control. Key: green nodes (circles) = model inputs; boxes = system state variables, thick arrows = flows, thin arrows = functional dependencies. All graphical notation is standard for ITHINK 7.0” software.
It collapses the multiple sub-compartments in the original model (needed to model age-structured cell populations and detailed maturation kinetics of proliferating and non-proliferating sub-populations) to only two aggregate compartments: a Reservoir of mature (non-proliferating) cells that are released to blood as needed; and a preceding compartment of Maturing proliferating cells that undergo successive rounds of cell division (amplification) and differentiation as they mature. This Maturing compartment, in turn, is supplied by stem cells, drawn from a deeper compartment that is not shown explicitly, since only the flux of new stem cells from it into the maturation process matters for modeling replenishment of the mature cell population.
The equations of the model, with their interpretations, are described next.
Reservoir of mature cells compartment
Interpretation
The mature cell reservoir provides a buffer of mature cells that are held in reserve to replace losses from the mature cell population as they occur. Such losses occur at a natural loss rate (indicated by the “normal_loss_rate” parameter), even in the absence of cytotoxic stress, as mature cells age and die. Cytotoxic exposures can increase the rate of losses from the mature cell compartment. Consistent with linear dynamics, the total loss rate, i.e., the “depletion” term, is modeled as directly proportional to the number of cells in the reservoir of mature cells. The constant of proportionality, i.e., the fractional loss rate per cell per unit time, is given by the sum: (normal_loss_rate + effective_dose), where the biologically “effective_dose” is defined and scaled (without loss of generality) so that an exposure producing an effective dose of 1 doubles the normal loss rate.
The reservoir of mature cells may be viewed as an inventory of replacement parts (cells) held to meet uncertain demands arising from unpredictable depletions of the mature cell population caused by environmental stresses. The feedback control system in Figure 1 represents the process that the tissue uses to try to maintain sufficient inventory to meet demands at all times. Thus, the feedback control law represents the tissue's strategy for managing the inventory of mature cells to protect against the risks of uncertain demands. The abstract problem of inventory control under uncertainty translates in physiological systems to the problem of homeostatic control of cell populations with kinetics linked via regulatory (e.g., cytokine signaling) networks. For example, for blood cells, the “Reservoir of mature cells” compartment refers specifically to a bone marrow reserve of mature granulocyte-macrophage cells that can be released to circulating blood as needed and that are replenished by proliferating cell populations within the marrow (Steinbach et al., 1980; Tibken and Hofer, 1995; Cox, 2000).
The initial size of the Reservoir of mature cells compartment can be normalized without loss of generality to have a value of 100 percent in unstressed steady-state equilibrium. The model will then describe the changing sizes of this cell population compartment in the presence of cytotoxic stresses as a percent of its normal (unstressed) value (i.e., 100% = normal value.) Similarly, the initial value of the Maturing compartment may also be normalized (without loss of generality) to 100% in unstressed steady-state equilibrium. Although multiple rounds of cell division (“amplification”) make the absolute number of mature cells exiting the differentiation/maturation pipeline several orders of magnitude greater than the number of proliferative cells entering it, scaling the compartments to express their sizes in terms of percentages of their respective normal equilibrium values (thus non-dimensionalizing the state variables) makes it unnecessary to estimate the absolute numbers of cells in the different populations for purposes of predicting relative (percentage) changes in response to cytotoxic exposures.
To reflect the fact that, in unstressed (zero-exposure) steady-state equilibrium, the rate at which replacement cells enter the Reservoir from the Maturing compartment must balance the routine loss of cells from the Reservoir due to normal depletion, the two numerical input constants maturation_rate and normal_loss_rate must be equal (since the normal equilibrium sizes of the Maturation and Reservoir compartments have both been normalized to 100%.) The common value of these two input parameters can be estimated empirically from average turn-over rates or cell residence times in either compartment. Differing numerical estimates from different cell kinetics experiments can potentially be reconciled and combined through generalized least squares or by Bayesian methods. For purposes of the conceptual illustrations in this paper, however, we simply assume that cells remain in the Reservoir for an average of about 20 days, so that maturation_rate = normal_loss_rate = 1/20 = 0.05 per day. For tissues with much slower or much faster kinetics than this roughly three-week average mature cell life implies, the time axes in our diagrams can be re-scaled to appropriately slow or speed responses, respectively. Thus, choosing a specific numerical value for these two parameters does not create any loss of generality (apart from a scaling constant for the time axis) in understanding the possible dynamic behaviors of the system. (Of course, as previously mentioned, the aggregation of all cell populations into only the few compartments shown in Figure 1 does deliberately over-simplify the detailed structure of the system, making it impossible to faithfully simulate detailed transients for arbitrary exposure inputs on a time scale of hours, as in Cox, 2000. But, as shown below, the simplified approximate model provides sharp insights into how non-linear responses can emerge even in linear dynamic models with delayed feedback control over cell populations.)
Maturing cells compartment
The kinetics of the Maturing (i.e., proliferating and differentiating) cells compartment are determined by the inflow of new stem cells into the maturation pipeline and the outflow from the pipeline into the Reservoir of mature cells. The maturation pipeline itself could be better modeled as a chain of sub-compartments (Cox, 2000) in order to more accurately describe the mean and variance of passage times through it. However, the key point for purposes of this paper is only that maturation creates an unavoidable delay between the need (“orders”) for replacement parts (cells) by the Reservoir of mature cells when its inventory becomes depleted and fulfillment of the need some time later (after stem cells recruited in response to the depletion of the Reservoir have had a chance to undergo amplification and differentiation into mature cells.) Placing a single Maturation compartment between new stem cells and the Reservoir compartment enforces such a delay (with an average value of (1/maturation_rate) = 20 days in this model). Even though the variance in the passage times is not well modeled by a single compartment (it is too large), the simple one-compartment model of maturation delay suffices to reproduce key features of the response dynamics from more complex and realistic models, such as over-compensating proliferation of cells following cessation of exposure (ibid). These are the features of central interest for understanding response nonlinearities to exposures that change only relatively infrequently (e.g., on a time scale with a length many times the average delay). Our computational experiments will be restricted to such exposure scenarios.
The outflow of cells from the Maturing compartment has already been described: it is just maturation = maturation_rate × Maturing, i.e., a constant fraction (e.g., 5% per day) of the cells in the Maturing compartment exit into the Reservoir_of_mature cells compartment per unit time, in order to replenish its losses. The final component of the system is the new_stem_cell piece, which provides the inflow to the Maturing compartment. We assume that, in the absence of exposure, stem cells are recruited and enter the maturation process at a rate sufficient to maintain the Maturing and Reservoir of mature cell compartments at their nominal levels of 100%. This is accomplished via a feedback loop that sends signals from the Reservoir compartment to govern the inflow to the Maturing compartment. [This inflow to the Maturing compartment is called “new_stem_cells”, not because the stem cells themselves are newly born—for example, they may have been lying quiescent in a G0 compartment for some time, possibly years or decades for tissues such as lung—but because they have newly entered the maturation pipeline that will result some time later (exponentially distributed with a mean of 1/maturation_rate = 20 days), after successive rounds of amplification and differentiation, in new mature cells flowing into the Reservoir compartment.] In unstressed steady-state equilibrium, all flows are equal (since all compartment sizes are normalized to 100), so new_stem_cells = maturation = depletion = maturation_rate × 100 = depletion_rate × 100 = 5.
When cytotoxic exposures deplete the Reservoir compartment more than normally, additional stem cells are recruited to enter the maturation pipeline to help compensate. This homeostatic process is modeled as follows.
Feedback Control Law
This feedback control law approximates physiological reality by positing that new stem cells enter the proliferation-differentiation (i.e., maturation) pipeline at a rate proportional to the depletion of cells below a nominal target level. The target_level input parameter is enough greater than the normal (zero-exposure equilibrium value) of 100 for the Reservoir of mature cells so that there is normally a non-zero flow of new stem cells (5 per day when all compartment values are at their nominal values of 100% with maturation and depletion flows of 5 per day). The greater the shortfall of the current actual level of the Reservoir of mature cells below the target level, the greater the recruitment of new cells to help close the gap. (If there is a surplus of mature cells above the target level, then no new stem cells are recruited.) Given a value for the feedback constant K, the value of the target_level parameter is determined from it by the requirement that compartment sizes should have their nominal values (of 100%) in the absence of exposure. Solving the steady-state equilibrium equation new_stem_cells = K × (target_level — 100) = 5 for target_level yields: target_level = (5/K) + 100. Thus, there is essentially only one parameter, K, to be estimated from data.
The constant K determines the strength of the feedback signal and the speed with which the system returns to equilibrium following a depletion of the Reservoir compartment. Its approximate numerical value can be identified from clinical and experimental data by choosing a value that reflects realistic adjustment speeds and qualitative patterns for specific tissues as well as possible using this simple model. A value of K = 1/7 = 0.1429 per day is used for the computational examples in this paper. (In terms of inventory control, this may be interpreted as follows: the system orders enough new stem cells per unit time to close the existing gap between actual and target levels of the Reservoir of mature cells in a week if there is no further depletion and if this order rate is maintained.) If much higher values are used (e.g., an order of magnitude larger), then transient exposures elicit qualitatively unrealistic responses with excessive oscillations as the system “over-steers” in its attempt to compensate for transient depletions of the Reservoir. If much smaller values are used (e.g., an order of magnitude smaller), then transient exposures have unrealistically long-lasting effects compared to experimental and clinical data (e.g., Cox, 2000; Tibken and Hofer, 1995). Thus, these rough qualitative constraints put useful bounds on the plausible numerical value of K.
For purposes of determining the possible shapes of dose-response relations, there is an important sense in which the entire system in Figure 1 involves only one unknown parameter. The normal sizes of the Maturing and Reservoir compartments need not be known, as we can normalize them to 100% and work in non-dimensionalized units (i.e., percentage changes). The maturation_rate and normal_loss_rate parameters are then equal and, since both express fractional rates per unit time, they can be assigned any arbitrary numerical value if the time axis is scaled correspondingly to get physiologically realistic response times. K determines the value of the target_level parameter by the formula: target_level = (5/K) + 100. Thus, the only input quantity that needs to be determined empirically is the ratio K to the maturation_rate parameter (e.g., 0.1429/0.05 = 2.858, in our examples.) This single parameter determines the shapes of the cell population trajectories in response to different exposure patterns in the simple model in Figure 1.
From Cytotoxicity to Possible Increased Cancer Risk
This section suggests a possible connection between cancer risk and cytotoxicity that stimulates the recruitment of abnormal numbers of stem cells (or other early, possibly multipotent, proliferative cells) into active cycling. Experimental evidence (e.g., Steinbach et al., 1980; Tibken and Hofer, 1995) shows that there are tight bounds on the rates at which cells can proliferate (e.g., the fraction of hematopoietic stem cells and early progenitor cells that proliferate instead of differentiating does not stray outside a range from about 0.4 to about 0.6 even under extreme stresses, while passage through the mitotic cycle duration can only be hastened to a very limited extent.) Therefore, substantial increases in the flux of stem cells entering the maturation pipeline can only be accomplished by recruiting additional stem cells that would normally be quiescent (e.g., in G0) into active cycling. If this reduces the time available to check and repair DNA damage, and/or necessitates recruiting very early, highly proliferative cells from deep compartments, the result may be a heightened risk of selecting initiated cells that can potentially undergo additional transformations to become carcinogenic.
Rather than attempting to model any of the details of this process, we will simply examine the dose-response implications if new_stem_cells flows above a certain threshold rate per unit time correspond to increased numbers of initiated cells per unit time entering a standard model of carcinogenesis (e.g., Hazelton et al., 2005.) We consider only cytotoxic challenges that do not kill the exposed animal, so that the size of the Reservoir compartment is never driven to zero, but stays at a substantial fraction of its normal value. Thus, we assume that enough new stem cells can always be recruited to satisfy demand (i.e., to fill all “orders” for new cells) to compensate for losses due to cytotoxicity, but that sufficiently high recruitment rates (“emergency orders”) increase cancer risks. (If all new stem cells that are recruited into actively cycling are at increased risk of carcinogenic damage, then the threshold may be zero.)
Figure 2 shows the structure of a cancer risk model in which values of the new_stem_cell flow above a certain threshold value are considered “emergency” orders that have increased probabilities of leading to cancer (e.g., they might contribute a flow of cells with excessive rates of unrepaired mutations to the “Initiated” compartment of a standard two-stage clonal expansion model of carcinogenesis). Their cumulative total is tracked over the course of the simulation via the “Cumulative emergency orders” variable. (Of course, cells from the Maturing compartment might also experience cytotoxic depletion and/or contribute to the flow of initiated cells. However, including these flows only complicates the model without revealing new model behaviors; hence, they are excluded for simplicity.) The top section of Figure 2 represents the same model as in Figure 1, but with a thin arrow drawn from normal_loss_rate to maturation_rate to indicate that these must have the same value when compartment sizes are normalized to 100% for normal equilibrium, and another arrow from K to target_level to indicate that the former determines the latter.
Structure of a cancer model with dose-dependent cytotoxicity.
Model Equation Summary
In summary, the complete set of model equations is as follows:
Input constants: R(0) = M(0) = 100, r = 0.05, K = 0.1429, R* = 135, threshold = 8.
Here, we use the abbreviations: R(t) = reservoir of mature cells at time t, M(t) = maturing cells at time t, and d(t) = effective dose at time t.
Design of Computational Experiments
Figure 2 also includes a dose model that allows two dosing intervals, separated by a recovery interval of duration specified by a spacing parameter between them, and a follow-up period after dosing ceases. Beginning at a user-specified start_time, a constant concentration (or dose rate) is applied for some duration (e.g., 6 hours). During this interval, mature cells are killed at a rate proportional to the concentration, where the constant of proportionality is a first-order cell-killing potency factor. This potency factor can be set equal to 1 without loss of generality, thus fixing the scaling for the effective dose axis. There is then a recovery interval of length spacing before the second dosing interval begins. Following the second dosing interval, all dosing ceases and the system is allowed to recover. Thus, the total cumulative effective dose received is 2 × C × T, where C = concentration and T = duration. [Although physiologically-based pharmacokinetic (PBPK) models could perhaps be used to explicitly model the conversion of administered doses to time courses of effective doses in relevant target tissues, the simple pattern of two dosing intervals separated by a gap suffices to describe the approximate internal dose profile corresponding to this pattern of administered doses, for a wide range of PBPK models (Cox, 1995).] Since the system is allowed to recover fully following the second exposure, the entire experiment can be repeated as many times as desired, thus multiplying the cumulative number of “emergency orders” and initiated stem cells produced by each replicate. The interesting dynamics take place within a single replicate, however (i.e., from initial equilibrium through the first dose, gap, second dose, and recovery period leading back to the initial equilibrium conditions.) Understanding the simulated dynamics of cell population responses to one such experiment is our goal, since multiplying by the number of replicates to calculate the cumulative effects of repeated cycles is then straightforward. (Of course, repeated high doses can permanently impair the ability of some tissues, such as the stromal microenvironment of the bone marrow, to respond effectively to future stresses. In this paper, as mentioned previously, we only consider relatively modest doses that do not come close to destroying cell populations and that allow full recovery following cessation of dosing.)
RESULTS
This section reports the results of several computational experiments carried out using the dynamic model and dosing experiment design in Figure 2. The main goal is to understand how nonlinear responses can arise even in a linear dynamic system with delayed feedback control. The response studied is the number of new stem cells recruited into the maturation process over time, which presumably might affect carcinogenesis by modulating the number of “normal stem cells” available at any time to be transformed (e.g., via unrepaired mutations) to initiated stem cells that enter the multistage stochastic process that may lead eventually to tumorigenic cells. (Mature cells are differentiated and so are not at risk of carcinogenic transformation.) Throughout the computational experiments reported here, the Reservoir of mature cells remains less than 135 (i.e., R(t) < R* in all experiments).
Results for C × T Experiments: C × T Is Not An Adequate Predictor of Response
The first set of experiments investigated how well the area-under-curve (AUC) of the effective dose predicts the excess risk caused by cytotoxic stress (assumed to be proportional to the cumulative “emergency orders” of stem cells if the dose is purely cytotoxic and has no effect on the carcinogenic process parameters, such as initiation, promotion, progression, repair, or apoptosis transition rates, apart from changing the number of proliferating stem cells available to enter the process.) The design was as follows. For each of several different C and T pairs having the same value of C × T, we examined the number of cumulative emergency orders created over the course of the experiment (using the commercial continuous simulation package ITHINK 7.0” to numerically solve the model). Table 1 summarizes the results. Clearly, the same C × T value for dose can elicit very different stem cell responses, depending on how the cumulative dose is allocated over time. A short, high dose has a disproportionately large effect compared to a longer-duration dose with the same AUC.
Different Cumulative Emergency Orders for the Same C × T Value.
C = concentration
T = duration
C × T
Cumulative Emergency Orders
0.01
32 days
0.32
0.00
0.04
8
0.32
0.65
0.08
4
0.32
5.13
0.16
2
0.32
11.70
0.32
1
0.32
24.36
Note: Spacing = 60 days for each of these dosing scenarios
Figure 3 shows the common-sense explanation revealed by the simulated time courses of new stem cell fluxes: unless the concentration is high enough to kill many mature cells quickly, the demand for new stem cells never spikes above the threshold level for “emergency orders” (set at 1.6 times the normal production rate of new stem cells in this example.)
Time courses of new stem cells for two dose scenarios with the same C × T.
Results on Spacing Between Exposures
The second set of computational experiments explores the role of the spacing parameter giving the time between the two dosing intervals. Of course, the value of spacing has no effect on the total dose (C × T) received by the tissue over the duration of the experiment. But it has a significant—and non-monotonic—effect on the stem cell kinetics response. Figure 4 shows the effect of different amounts of spacing, from 5 days to 110 days, between successive doses with a concentration of C = 0.08 for a duration of T = 4 days.
The Time (“Spacing”) Between Doses Affects Stem Cell Responses.
Although increasing the recovery interval between doses is beneficial at first (up to about 25 days), further increases (out to about 60 days) actually increase response. Beyond about 90 days, the two doses are so widely separated in time that they effectively do not interact—the cumulative response of emergency orders is just twice the response to either of the two doses. Inspecting the detailed trajectories of the new stem cell fluxes into the “Maturing” compartment for different values of spacing shows that this surprising non-monotonic pattern arises from the crowding together in time of responses to the two dose intervals. If the second dose occurs when the system has only partly responded to the first, then the effect of the second dose may be either greater or lesser than its effect after full adjustment has taken place, depending on how much of the response to the first dose has occurred when the second dose begins. Figures 3 and 4 each imply that smaller total doses may have larger cytotoxic effects on stem cells, depending on how the doses are allocated over time. However, the explanatory mechanisms involved differ: nonlinear dependence of peak stem cell production on C × T in Figure 3vs. nonlinear dependence of cumulative peak stem cell production on the duration of the zero-exposure interval between two positive dosing intervals in Figure 4.
Results for a Model with Exposure-Induced DNA Repair
Rather than only counting the total number of emergency orders, it is useful to consider the fraction of such stem cells that are likely to become initiated cells due to unrepaired DNA damage. To model the fact that some initiated cells occur even in the absence of exposure, the threshold parameter in the model can be reduced to less than the normal level of new_stem_cells, so that a small flow of emergency orders occurs even without exposure. Exposures that induce protective enzymes or other mechanisms that prevent or repair initiation damage in cells before it is locked in by mitosis may then reduce risk. For example, suppose that the background initiation rate is A per “emergency order” cell per unit time and that the repair rate is (B + Q × effective_dose) per initiated cell per unit time, where B is a background (dose-independent) repair rate and Q is a dose-induced repair potency factor. Then the probability that mitosis locks in unrepaired damage, denoted by p, can be approximately modeled as the outcome of an equilibrium between the rate of initiation and the rate of repair, satisfying: A(1 — p) = (B + Q × effective_dose) p, i.e., in equilibrium, the probability that a cell has unrepaired damage times the repair rate equals the probability that it is undamaged times the damage rate. (The competing initiation and repair processes can be represented in more detail via an additional pair of ordinary differential equations, but if the time scale for repair is fast compared to the time scale for cell population kinetics, then the equilibrium can be used as an approximation to the explicit dynamic process.) This yields the formula p = A/(A + B + Q × effective_dose), which is equivalent to: p = 1/(b + q × effective_dose), where b and q are “reduced parameters” that might be estimated from data (or from the identities b = (1 + B/A) and q = Q/A, if the “structural parameters” A, B, and Q are known.) If repair is highly reliable even without exposure, then b will be large and p will be small (compared to 1). Induced repair or other mechanisms that reduce p in the presence of exposure cannot reduce it below zero, no matter how large q is. Thus, the fraction of unrepaired “emergency order” stem cells will be at most p = 1/b (while there is no exposure) and at least p = 0 (even during exposure). (For large q, this step function provides an excellent approximation to the function p.)
Figure 5 shows the result of a simulation run in which the probability of unrepaired damage when an “emergency order” stem cell undergoes mitosis is given by: p = 1/(100 + 100000 × effective_dose), i.e., for b = 100, q = 100,000. [The y axis is scaled so that “1” represents that number of initiated cells (unrepaired “emergency orders”) formed in one year when the value of the threshold parameter is 4 (and duration = 7, spacing = 20. The absolute value of this number is 3.65 initiated cells per year.)] The initial 2.5% reduction in risk is caused by the fact that even low concentrations induce a protective effect (e.g., metallothionein induced by cadmium.) At higher concentrations, the stem cell proliferation induced by cytotoxicity leads to elevated “emergency order” rates that persist after exposure ceases, leading to a net increase in risk. This interaction between induced protection and stem cell proliferation as tissue populations return to normal provides a possible explanation for U-shaped dose-response relations different from the better-known apoptosis-mediated mechanism, in which exposure-related cell-killing of initiated populations provides a net protective effect over some range of doses (Bogen, 1998; Holt, 1997 for radiation.)
Induced Repair of Initiation Damage Can Produce Hormesis (U-Shape) [duration = 7 days, spacing = 20 days, p = 1/(100 + 100000 × effective_dose)].
Figure 6 shows the results of re-running the simulation in Figure 6 with the following changes in inputs: the duration of each dose is extended from 7 days to 40 days and p = 1/(1 + 10 × effective_dose), i.e., damage repair is only induced only when exposure is present. In this case, rather than being U-shaped, the dose-response relation exhibits an n-shape, which might be termed inverse hormesis: only when concentrations are sufficiently large (and exposure durations are sufficiently long) does a protective effect occur due to induced repair outweighing stem cell proliferation.
This paper has examined dose-response relations for cytotoxicity and cytotoxicity-mediated carcinogenesis based on modulation of the number of stem cells entering the carcinogenic process. In contrast to several earlier papers on hormesis in carcinogenesis models (e.g., Bogen, 1998), we treat the entire carcinogenic process as a black box and focus not on what might happen inside it—such as effects of exposure on proliferation and apoptosis of initiated cells (promotion) or transformation of initiated to malignant cells (progression)—but rather on how many cells per unit time are expected to enter the process.
Our major conclusions are that: (a) Even when normal stem cells dynamics are described by an underlying linear dynamic system (represented by a system of constant-coefficient linear ordinary differential equations) with linear feedback control (new stem cells enter the differentiation/maturation pipeline at a rate proportional to the gap between target and actual levels, which is always positive in our analyses), still the quantities of greatest interest for risk assessment, namely, “emergency orders” of new stem cells produced per year, can be strongly nonlinear (and non-monotonic) functions of dose rate. For example, in the system shown in Figures 1 and 2, even continuous exposures to constant concentrations less than 0.02 will not create any emergency orders. (b) The total dose, C × T, received over a given interval does not provide a sufficient basis for predicting (or even approximating) resulting risk (see Table 1 and Figure 4.) Instead, how the cumulative dose is allocated over time can greatly affect response. (c) If low levels of exposure induce a protective response that reduces the fraction of stressed-production stem cells (“emergency orders”) that become initiated cells entering the carcinogenesis sub-model, then a U-shaped or J-shaped dose-response curve (Figure 5) or an n-shaped dose-response curve (Figure 6) can result. For the J-shaped curve in Figure 5, higher concentrations eventually produce enough emergency orders per unit time to overwhelm the protective effects achieved at lower concentrations. By contrast, the n-shaped dose-response relation in Figure 6 results when exposure duration is made long enough so that the protective effect induced by exposure outweighs the additional emergency orders created as the system compensates for depletion of mature cells by cytotoxicity. (As previously mentioned, proliferating cells could also be killed by cytotoxic exposures, slightly complicating the model, but simulations that include this added complexity produce qualitative results similar to those already illustrated.)
The cell inventory-management perspective illustrated in this paper is designed to complement other models of mechanisms of hormesis for two-stage clonal expansion (TSCE) models of carcinogenesis (e.g., Hazelton et al., 2005; Bogen, 1998; Holt, 1997). Thus, we have focused on carcinogenic effects that are not due to familiar mechanisms of direct genotoxic initiation, promotion, or progression. None of the standard parameters of a TSCE model, such as the initiating mutation rate, the net birth rate of initiated cells, or the rate of transformation of initiated to malignant cells, need be affected by dose to produce the results we have discussed. Rather, we have considered only the possible dose-response effects due to compensating proliferation of stem cells (and, in Figures 5 and 6, exposure-induced stimulation of damage repair) in response to cytotoxic exposures. Some important chemicals appear to be carcinogenic only when they induce compensating proliferation of normal cells (Cox, 2005), so this modeling focus may be of practical interest. A promising future research direction is to combine such modeling of purely cytotoxic effects on normal stem cells and cancer risks with previous models of effects of carcinogenic exposures on TSCE model parameters.
Finally, although this paper has only considered applications to cell populations, any system that satisfies similar equations describing homeostatic control of a reservoir of a finished resource (e.g., protective enzymes, intracellular proteins, etc.) held to meet sudden, unpredictable demands and maintained via delayed feedback control over production can be expected to exhibit nonlinearities similar to those described here.
Footnotes
ACKNOWLEDGEMENTS
The author gratefully acknowledges a research grant from the American Chemistry Council in 2002–2003 (American Chemistry Council Project: RSK-01-03) that supported my initial research on simulation modeling of cytotoxic effects on stem cell compartments;and support from Philip Morris International in 2004–2005 for my research on computational models of lung carcinogenesis,which allowed development of the framework and results presented here.
References
1.
AtkinsHL. A clinician's comments on the cyclophosphamide hematotoxicity biologically based risk assessment model. J. Toxicol. Environ. Health A.2000 Nov;61(5–6): 529–38.
2.
BogenKT. Mechanistic model predicts a U-shaped relation of radon exposure to lung cancer risk reflected in combined occupational and US residential data. Hum. Exp. Toxicol.1998 Dec;17(12): 691–6.
3.
CalabreseEJ. Toxicological awakenings: The rebirth of hormesis as a central pillar of toxicology. Toxicol. Appl. Pharmacol.2005 Apr 1;204(1): 1–8.
4.
CoxLAJr. Universality of J-shaped and U-shaped dose-response relations as emergent properties of stochastic transition systems. Nonlinearity in Biology, Toxicology, and Medicine3: 285–198. 2005.
5.
CoxLAJr. A biomathematical model of cyclophosphamide hematotoxicity. J. Toxicol. Environ. Health A.2000 Nov;61(5–6): 501–10.
6.
CoxLAJr. Simple relations between administered and internal doses in compartmental flow models. Risk Anal.1995 Apr;15(2): 197–204.
7.
FliednerTMTibkenBHoferEPPaulW. Stem cell responses after radiation exposure: A key to the evaluation and prediction of its effects. Health Phys.1996 Jun;70(6): 787–97.
8.
FukushimaSKinoshitaAPuatanachokchaiRKushidaMWanibuchiHMorimuraK. Hormesis and dose response-mediated mechanisms in carcinogenesis: Evidence for a threshold in carcinogenicity of non-genotoxic carcinogens. Carcinogenesis2005 Jun 23
9.
HazeltonWDClementsMSMoolgavkarSH. Multistage carcinogenesis and lung cancer mortality in three cohorts. Cancer Epidemiol. Biomarkers Prev.2005 May;14(5): 1171–81.
10.
SteinbachKHRafflerHPabstGFliednerTM. A mathematical model of canine granulocytopoiesis. J. Math. Biol.1980 Aug;10(1): 1–12.
11.
TibkenBHoferEP. A biomathematical model of granulocytopoiesis for estimation of stem cell numbers. Stem Cells1995 May;13Suppl 1: 283–9.
