CarcinoMod: a tool for stochastic carcinogenesis modeling

Qi Zheng
National Center for Toxicological Research
Jefferson, Arkansas 72079

Introduction

Since as early as 1960,  stochastic birth-and-death processes have been the predominant building blocks for constructing biologically based carcinogenesis models (Kendall 1960). The two-stage carcinogenesis theory was then succinctly  couched in the birth-and-death process language by Kendall's pioneering work, resulting in three types of stochastic carcinogenesis models. Moolgavkar and Venzon's (1979) subsequent seminal work recast the two-stage models as important tools for connecting  cancer survival and hazard rate functions with the biologically  most important cellular kinetic parameters, i.e., cell division rate, cell death rate and cell mutation rate. The survival and hazard rate functions are sometimes justifiably called two fundamental quantities, as these two quantities were  inextricably woven into the majority of carcinogenesis modeling efforts during the last two decades or so. 

However, compelling biological evidence suggests that two mutations are not sufficient to cause certain cancers. For example,  Kopp-Schneider and Portier (1995) analyzed mouse skin painting experiment data that indicate the need for more than two stages in cancer models. Regarding colorectal tumor, Fearon and Vogelstein (1990) commented: "mutation in at least four or five genes are required for the formation of a malignant tumor." Using radiation effects data, Little (1995) also argued that some extra stages should be appended to the two-stage cancer models. At the instigation of those authors,   Zheng (1997a) offered systematic computational recipes for a rather broad class of multistage carcinogenesis models. The package CarcinoMod is a Mathematica implementation of all those recipes. Specifically, CarcinoMod provides functions for computing the survival and hazard rate functions for all eight types of models presented in Zheng (1997a). CarcinoMod can handle  models with both constant parameters and time-dependent parameters. Moreover, for the extensively studied (Subscript[B, 2])/0/0 model, all existing closed form expressions for the survival and hazard rate functions are separately implemented, enhancing computational accuracy and efficiency for this particularly simple model.

The goal of this notebook is to demonstrate the various functionalities that CarcinoMod provides. For  mathematical details, the reader is referred to the technical papers listed in the reference section. We begin with the simpler cases --- models with constant parameters.

SetDirectory["~/CarcinoMod"]

<< CarcinoMod.m

General k-stage model with constant parameters

For models with constant parameters (also called homogeneous models), the name of a CarcinoMod function  is composed of two parts: the model name and the quantity to be computed. For example, to compute the survival function for a (Subscript[B, k])/N/0 model, we invoke the function BN0Survival. The first part of the function name "BN0" indicates that this function is for (Subscript[B, k])/N/0 models where k is an arbitrary positive integer. The function name itself does not contain any information about the exact number of stages in the model, because k can be inferred by checking the number of cellular parameters supplied by the user. The second part "Survival" is simply an English word with its first letter capitalize, the meaning of which is self-explanatory. By analogy, the function name A00Hazard suggests that this function is for computing the hazard rate function for  (Subscript[A, k])/0/0 models. It is particularly worth noting that, for any constant parameter model considered in Zheng (1997a), the solution of one ordinary differential equation (ODE) system yields all the values for both the survival and the hazard rate functions. Therefore, these functions are so designed that, by specifying a positive real number T, the user get  the survival (or hazard rate) function as an interpolating function object that are valid on the whole interval [0,T].

?BN0Survival

The following code draws the survival function for a (Subscript[B, 3])/N/0 model on the time interval [0,2000]. The process starts with five normal cells.

surv = 
      BN0Survival[
        2000, {{0.07, 0.063, 0.0001}, {0.068, 0.062, 0.0001}, {0.07, 0.063, 
            0.0001}}, 5]

Plot[surv[t], {t, 0, 2000}, Frame -> True]

Objects returned by CarcinoMod as interpolating function objects  can just be treated as  "ordinary" functions, as long as the argument is within its  valid domain . For instance, the survival probability at time t=1234 for the above model was already "stored" in the interpolating function returned by BN0Survival.

surv[1234]

We now plot the hazard rate function for an (Subscript[A, 3])/N/0 model, using the same  parameter values.

haz = 
      AN0Hazard[
        2500, {{0.07, 0.063, 0.0001}, {0.068, 0.062, 0.0001}, {0.07, 0.063, 
            0.0001}}, 5]

Plot[haz[t], {t, 0, 2500}, Frame -> True]

Before introducing more advanced examples, we make an informal assessment of the general computational approach of CarcinoMod as described in Zheng (1997a). Because closed -form formulae for the (Subscript[B, 2])/0/0 model are available  (Kopp-Schneider et al. 1994 and Zheng 1994) and was separately built into  CarcinoMod, we can use the (Subscript[B, 2])/0/0 model as a benchmark.

?B200Survival

survExact[t_] := 
      B200Survival[t, 5.5, {23.7,  23.58, 5.5*10^(-7)}]

g1 = 
      Plot[survExact[t], {t, 0, 250}, PlotStyle -> RGBColor[0, 1, 0], 
        DisplayFunction -> Identity];

survODE = 
      B00Survival[250, 5.5, {{23.7, 23.58, 5.5*10^(-7)}}]

g2 = 
      Plot[survODE[t], {t, 0, 250}, PlotStyle -> RGBColor[1, 0, 1], 
        DisplayFunction -> Identity];

Show[g1, g2, Frame -> True, DisplayFunction -> $DisplayFunction]

The two survival curves produced using the same (Subscript[B, 2])/0/0 model but with different computational approaches are almost indistinguishable. It indicates that the ODE based general approach performs satisfactorily, at least in this example. The two corresponding hazard rate curves also match excellently.

hazExact[t_] := 
      B200Hazard[t, 5.5, {23.7,  23.58, 5.5*10^(-7)}]

g1 = 
      Plot[hazExact[t], {t, 0, 120}, PlotStyle -> RGBColor[0, 1, 0], 
        DisplayFunction -> Identity];

hazODE = B00Hazard[120, 5.5, {{23.7, 23.58, 5.5*10^(-7)}}]

g2 = 
      Plot[hazODE[t], {t, 0, 120}, PlotStyle -> RGBColor[1, 0, 1], 
        DisplayFunction -> Identity];

Show[g1, g2, Frame -> True, DisplayFunction -> $DisplayFunction]

Models with time-dependent parameters

Using models with time-dependent parameters (also called nonhomogeneous models) in general requires considerably more computational resources. Specifically, for each given (Subscript[t, 0]) > 0, one needs to solve an ODE system to obtain the survival probability S((Subscript[t, 0])). For this reason, each CarcinoMod function for computing survival probabilities for nonhomogeneous models returns only numerical quantities, instead of interpolating  function objects. Time dependency of a particular parameter is specified by a pure function. For example, if a cell birth rate varies with the time t as β(t) = 0.028 (1 + 0.00001 (t^2)),  it is supplied to a CarcinoMod function either as 0.028 (1+0.00001 #)& or as Function[t, 0.028 (1+0.00001 t^2)]. The dummy variable t can be replaces with any  legitimate symbols. The name of a CarcinoMod function for a nonhomogeneous model has the same structure as the name for its corresponding homogeneous model, except that the second part of the name always begins with "Nonhomo". For example, the function BN0NonhomoSurvival is the counterpart function of BN0Survival.

?BN0NonhomoSurvival

Consider a (Subscript[B, 2])/N/0 model. Assume that the cellular birth, death and mutation rates for normal (stage 0) cells are \(TraditionalForm
      \`\(β\_0\)(t) = 0.03\ \((1 + 0.00000001\ t\^2)\)\), (Subscript[δ, 0])(t) = 0.028 (1 + 0.00001 t)  and (Subscript[μ, 0])(t) = 0.00001, and for initiated (stage 1) cells  are (Subscript[β, 1])(t) = 0.03, (Subscript[δ, 1])(t)=0.028, and (Subscript[μ, 1])(t) = 0.00001. The survival probability at t=12345 can be found by calling BN0NonhomoSurvival.

BN0NonhomoSurvival[
      12345, {{0.030 (1 + 0.00000001  #^2)&, 
          0.028 (1 + 0.00001 #)&, 0.00001&}, {0.030&, 0.028&, 
          0.00001&}}, 5]

Although pure functions provide a very flexible way of supplying time-dependent parameters to CarcinoMod functions, it is sometimes desirable that parameters be given in traditional form. For example, if all but only one or two parameters are constant quantities, it is tedious to type a long list of parameters as pure functions. The function MakeParameterFunctional allows the user to type parameters in traditional form and then transform them into pure function.

?MakeParameterFunctional

Clear[u]; 
    para = {{0.028 (1 + Sin[0.01  u]),  0.025, 0.0001}, {0.028, 0.020, 
          0.002}}

para0 = MakeParameterFunctional[para, u]

surv = A00NonhomoSurvival[#, 0.02&, para0]&

Plot[surv[t], {t, 0, 1234}, Frame -> True]

We  recall briefly that, if S(t) denotes a survival function, the corresponding hazard rate function is defined by  h(t) = -S' (t)/S(t) where S' (t) is the derivative of the survival function S(t). Because the hazard rate  at a given point (Subscript[t, 0]) ≥ 0 reflects the behavior of the corresponding survival function in a neighborhood of that point, it is necessary to compute a number of survival probabilities around (Subscript[t, 0]) in order to obtain the hazard rate at (Subscript[t, 0]). Consequently, with nonhomogeneous models, several ODE systems must be solved to get one value of the hazard rate function. It is therefore natural and economical that, with time-dependent parameter models, CarcinoMod functions return interpolating function objects for hazard rate functions, as in the case of homogeneous models. Numerical differentiation of the survival function S(t) is simple in principle but tricky in practice, if the survival function is not sufficiently smooth. In order to compute S' (t), CarcinoMod first samples the survival function S(t), creates an interpolation function object for S(t), and then differentiates the interpolating function. Three options are available that allow the user to control how these interpolating function objects are constructed. The option HazardInterpolationPoints specifies how many sampling points should be used to compute the corresponding survival probabilities. The second option, HazardInterpolationOrder, specifies the order of the interpolating function that will eventually be returned to the user. The third option, HazardBoundaryStretch, if set to some ϵ > 0, will enlarge the prescribed interval from [0,T] to [0, (1 + ϵ) T]. The purpose of extending the prescribed time interval [0,T] is to reduce the boundary effect of interpolation and hence enhance accuracy around the boundary point T.

HazardOptions[AN0NonhomoHazard]

SetHazardOptions[AN0NonhomoHazard, HazardInterpolationPoints -> 50]

haz = 
      AN0NonhomoHazard[
        3000, {{0.028 (1 + 0.0002 #)&, 0.025 &, 0.00001&}, {0.028&, 
            0.025&, 0.00001&}}, 10, HazardInterpolationOrder -> 16, 
        HazardInterpolationPoints -> 100, HazardBoundaryStretch -> 0.1]

haz[1500]

Plot[haz[s], {s, 0, 3000}, PlotRange -> All, Frame -> True]

Functions for the (Subscript[B, 2])/0/0 model

Due to its structural simplicity, the (Subscript[B, 2])/0/0 model possesses several closed-from analytic formulae for the survival and hazard rate functions that were published during the past few years. In fact, all the analytic formulae derived for the survival and hazard rate functions for the (Subscript[B, 2])/0/0 models are also valid for the (Subscript[A, 2])/0/0 and (Subscript[C, 2])/0/0 models (see Zheng 1998). To retain the historical flavors of the field, CarcinoMod functions  using these analytic formulae all begin with the prefix "B200" as if they were solely for the (Subscript[B, 2])/0/0 models. We have already encountered two of these functions: B200Survival and B200Hazard. Another function in this group is  B200PiecewiseLinearHazard which computes the hazard rate function for a (Subscript[B, 2])/0/0 model that allows the first mutation rate (the intensity of a Poisson process) to change with the time in a piecewise linear manner. The formula was derived by this author in an unpublished manuscript in 1995 and is useful when the tissue under consideration is still growing. Because in a (Subscript[B, 2])/0/0 model,  normal cell growth is assumed to follow some deterministic law (e.g., the logistic law) and  the rate at which normal cells mutate to initiated cells is proportional to the size of the normal cell population, it is  not reasonable to assume constant initiation rates when dealing with  growing tissues. Piecewise linear functions allow the user to approximate the normal cell  growth pattern with piecewise linear functions.

?B200PiecewiseLinearHazard

In the following example, we assume that a growth function is composed of three pieces of linear functions.

normalGrowth = {{200, 0.01, 0.00195}, {400, 0.25, 0.00075}, {Infinity, 
          0.51, 0.0001}}

For comparison purposes, we first draw the graph of this piecewise linear function. We will assume in the subsequent example that the total number of  normal cells susceptible to mutation at  any given time t is proportional to the value of this piecewise linear function at t.

ListPlot[{{0, 0.01}, {200, 0.4}, {400, 0.55}, {1000, 0.61}}, 
      PlotJoined -> True, Frame -> True]

Now, we further assume that the birth, death and mutation rates for initiated cells are 0.028, 0.027 and 0.0001, respectively. The hazard rate function based on a (Subscript[B, 2])/0/0 model can now be easily drawn by invoking B200PiecewiseLinearHazard.

haz = 
      B200PiecewiseLinearHazard[#, normalGrowth, {0.028, 0.027, 0.0001}]&

Plot[haz[t], {t, 0, 1000}, Frame -> True]

The last two functions in this group deal with the (Subscript[B, 2])/0/0 model with all four parameters piecewise constant. The algorithm for computing the survival function for this case can be found in Zheng (1995a). It was first implemented with Mathematica as a user function for MathSource item  0207-931 (Zheng 1995b).   The algorithm for computing the hazard rate function for the same scenario was recently published by Heidenreich et al. (1997), and is also a user function of CarcinoMod. We here present the same five-piece example from MathSource item 0207-931.

\(\(fivePieces = {{4000, 0.0017, 0.0015, 0.03, 0.0000003}, 

        \ \ \ \ \ \ \ \ \ \ \ \ {8000, 0.0015, 0.0014, 0.02, 0.0000001}, 

        \ \ \ \ \ \ \ \ \ \ \ {12000, 0.0016, 0.0016, 0.03, 0.0000002}, 

        \ \ \ \ \ \ \ \ \ \ \ {16000, 0.0014, 0.0012, 0.02, 0.0000005}, 

        \ \ \ \ \ \ \ \ {Infinity, 0.0015, 0.0013, 0.01, 0.0000004}}; \)\)

?B200PiecewiseConstantSurvival

surv = B200PiecewiseConstantSurvival[#, fivePieces]&

Plot[surv[t], {t, 0, 20000}, Frame -> True, GridLines -> Automatic]

haz = B200PiecewiseConstantHazard[#, fivePieces]&

Plot[haz[t], {t, 0, 20000}, Frame -> True, GridLines -> Automatic]

Concerning the last two functions, we  offer a caveat. Even if the exposure dose is perfectly piecewise constant, the internal dose is almost always continuous, as already demonstrated by many physiologically based pharmacokinetic (PBPK) models. See Zheng (1997b) for a PBPK model implemented with Mathematica. Caution therefore must be used before adopting carcinogenesis models with piecewise constant parameters.

Summary

CarcinoMod provides the user with 32 functions to compute the survival and hazard rate functions for the eight types of carcinogenesis models given in Zheng (1997a). For each type of models, there are four user functions whose names have a common 3-letter prefix to indicate the models that the functions address. For example, the prefix "BN0" signifies functions for the (Subscript[B, k])/N/0 models. This prefix is followed by either Survival, NonhomoSurvival, Hazard or NonhomoHazard. Such function names are self-explanatory and easy to remember. All these 32 functions adopt the ODE-based computational approach as described in Zheng (1997a). In addition, there are five functions whose names have the prefix "B200". These five function names no doubt  remind the user of the (Subscript[B, 2])/0/0 model, although they are in fact valid for the (Subscript[A, 2])/0/0 and (Subscript[C, 2])/0/0 models as well. These five functions use exact formulae to compute the survival and hazard rate functions. Many of the stochastic carcinogenesis models developed during the last two decades or so are sequential, unidirectional compartmental models, and all such models can be handled by CarcinoMod.

Acknowledgments

A major part of this work was performed while the author was at Wolfram Research, Inc. as a Mathematica visiting scholar in December, 1997.

References

Fearon, E.R. and Vogelstein, B. (1990). A genetic model for colorectal tumorigenesis, Cell, 61:759-767.

Heidenreich, W.F., Luebeck, E.G. and Moolgavkar, S.H. (1997). Some properties of the hazard function of the two-mutation clonal expansion model, Risk Analysis, 17:391-399.

Kendall, D.G. (1960). Birth-and-death processes, and the theory of carcinogenesis, Biometrika, 47:13-21.

Kopp-Schneider, A., Portier, C.J. and Sherman, C. D. (1994). The exact formula for tumor incidence in the two-stage model, Risk Analysis, 14:1079-1080.

Kopp-Schneider, A. and Portier, C.J. (1995). Carcinoma formation in NMRI mouse skin painting studies is a process suggesting greater than two stages, Carcinogenesis, 16:53-59.

Little, M.P. (1995). Are two mutations sufficient to cause cancer? Some generalizations of the two-mutation model of carcinogenesis of Moolgavkar, Venzon, and Knudson, and of the multistage model of Armitage and Doll, Biometrics, 51:1278-1291.

Moolgavkar, S.H. and Venzon, D.J. (1979). Two-event models for carcinogenesis: Incidence curves for childhood and adult tumors, Mathematical Biosciences, 47:55-77.

Zheng, Q. (1994). On the exact hazard and survival functions of the MVK stochastic carcinogenesis models, Risk Analysis, 14:1081-1084.

Zheng, Q. (1995a). On the MVK stochastic carcinogenesis model with Erlang distributed cell life lengths, Risk Analysis, 15:495-502.

Zheng, Q. (1995b). MathSource item 0207-931, abstracted in The Mathematica Journal Volume 6 (1996), Issue 1, p.17 .

Zheng, Q. (1997a). A unified approach to a class of stochastic carcinogenesis models, Risk Analysis, 17:617-624.

Zheng, Q. (1997b). Exploring physiologically based pharmacokinetic models, Mathematica in Education and Research, 6:22-28.

Zheng, Q. (1998). On a trinity role of the hazard and survival functions of the two-stage carcinogenesis models, Communications in statistics -- Computation and Simulation, to appear.