CarcinoMod: A Tool for Stochastic Carcinogenesis Modeling

Qi Zheng
National Center for Toxicological Research
Jefferson, Arkansas 72079

Introduction

Stochastic birth-and-death processes have been the predominant building blocks for constructing biologically based carcinogenesis models for several decades (Kendall 1960). The two-stage carcinogenesis theory was  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 functions with  biologically  important cellular kinetic parameters, i.e., cell division rate, cell death rate and cell mutation rate. The survival and hazard 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. 

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 indicated 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. Inspired by these developments,  Zheng (1997a) refined and generalized an algorithm  proposed by Quinn (1989) for computing survival functions for the (Subscript[B, 2])/0/0 model and  offered systematic computational recipes for a rather broad class of multistage carcinogenesis models. The package CarcinoMod is a Mathematica implementation of all those recipes.  CarcinoMod provides the user with a complete set of  functions to compute the survival functions, the  hazard functions and the cumulant functions for the eight types of carcinogenesis models given in Zheng (1997a). For each type of models, there are six user functions whose names have a common 3-letter prefix to indicate the type of 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 one of the following six self-explanatory shorthand names:  Survival, NonhomoSurvival, Hazard, NonhomoHazard, Cumulant and CumulantODE. These function names are  easy to remember. For instance, the function BN0Cumulant can be invoked to compute a cumulant function for the (Subscript[B, k])/N/0 model.  To enhance computational accuracy and efficiency for the extensively studied (Subscript[B, 2])/0/0 model, all existing closed form expressions for the survival and hazard rate functions for this model were separately implemented, resulting in five functions whose names have the prefix "B200". In addition,  CarcinoMod also contains 2 auxiliary user functions to facilitate the use of the package for certain special cases. All told, CarcinoMod offers  67 user functions.

The goal of the present notebook is to demonstrate the various functionalities that CarcinoMod provides. For an  overview of the field, the reader can consult Moolgavkar and Luebeck (1990), Tan (1991), and Kopp-Schneider (1997). 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 time-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 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 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

Denes and Krewski (1996) presented graphs for the survival and hazard functions for two (Subscript[B, 2])/N/0 models. We reproduce their Figure 1 for Example I (the exact solutions only). The following code draws the survival function.

surv = 
      BN0Survival[60, {{1, 0.6, 10^(-6)}, {1, 0.3, 10^(-6)}}]

Plot[surv[t], {t, 0, 60}, PlotRange -> All, Frame -> True]

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

surv[55]

The hazard function can be similarly drawn.

haz = 
      BN0Hazard[60, {{1, 0.6, 10^(-6)}, {1, 0.3, 10^(-6)}}]

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

To create a model of more stages, one only needs to supply  more cellular kinetic parameters. For example, the hazard function of  a (Subscript[B, 4])/N/0 model can be drawn as follows.

haz4stage = 
      BN0Hazard[
        120, {{1, 0.6, 10^(-6)}, {1, 0.6, 10^(-6)}, {1, 0.3, 
            10^(-6)}, {1, 0.3, 10^(-6)}}]

Plot[haz4stage[t], {t, 0, 120}, Frame -> True, PlotRange -> All]

Portier et al. (1997) compared the tumor probability functions between a (Subscript[B, 2])/N/0 model and a (Subscript[C, 2])/N/0 model. Their Figure 3 can be reproduced as follows.

sB = 
      BN0Survival[
        80, {{100, 100, 4.5*10^(-7)}, {103.5, 103.29, 
            4.5*10^(-7)}}, 10^8]

gB = 
      Plot[1 - sB[t], {t, 0, 80}, PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction -> Identity];

sC = 
      CN0Survival[
        80, {{100, 100, 4.5*10^(-7)}, {103.5, 103.29, 
            4.5*10^(-7)}}, 10^8]

gC = 
      Plot[1 - sC[t], {t, 0, 80}, PlotStyle -> RGBColor[0, 1, 0], 
        DisplayFunction -> Identity];

Show[gB, gC, Frame -> True, DisplayFunction -> $DisplayFunction];

We now 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 function curves produced using the same (Subscript[B, 2])/0/0 model but with different computational approaches are almost indistinguishable. The two corresponding hazard function 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  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)\)\), \(TraditionalForm
      \`\(δ\_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 specifying time-dependent parameters, 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];

A biologically  meaningful way of choosing time-dependent model parameters is to incorporate  physiologically based pharmacokinetic (PBPK) models into stochastic carcinogenesis models.  This approach was considered as "the most promising application of carcinogenesis models for risk assessment" by Kopp-Schneider (1997).  We give a simple example based on  Conolly et al. (1988) who combined the PBPK model of Ramsey and Andersen (1984) with a (Subscript[B, 2])/0/0 carcinogenesis model.  First we need produce the ODE system that describes kinetics of  styrene  inhaled by a laboratory rat (the code below is adapted from  Zheng 1997b).

PBPKode={AmtFat'[t]==0.391134  ConcInhale - 0.338178  AmtFat[t] + 
                5.660445  AmtLiver[t] + 0.269885  AmtMuscle[t] + 
                2.38899  AmtRichly[t],
		AmtMuscle'[t]==0.531943  ConcInhale + 0.0437813  AmtFat[t] + 
                7.69821  AmtLiver[t] - 2.7379798  AmtMuscle[t] + 
                3.24903  AmtRichly[t],
		AmtRichly'[t]==1.83833  ConcInhale + 0.1513023  AmtFat[t] + 
                26.6041  AmtLiver[t] + 1.26846  AmtMuscle[t] - 
                16.2571  AmtRichly[t],
		AmtLiver'[t]==1.65059  ConcInhale + 0.135851  AmtFat[t] + 
                (-41.23634  + 
                      111.111/(-0.36  - 30.8642  AmtLiver[t])) 
                  AmtLiver[t] + 1.13891  AmtMuscle[t] + 
                10.0816  AmtRichly[t]};
\(PBPKinit = {AmtFat[0] == 0, AmtMuscle[0] == 0, 
AmtRichly[0] == 0, 
          AmtLiver[0] == 0}\);

We then assume that the styrene exposure  can be described by a step function.

ppmToMgl[ppm_]:= ppm 104.151/10^3/0.0821/298;
ExposurePattern={ConcInhale->If[t<20.0||(t>40&&t<60),ppmToMgl[1200.],0.]};
FinalRatODE=Join[PBPKode,PBPKinit]/.ExposurePattern//Simplify;

The ODE system for  the PBPK model must be solved numerically before the PBPK model can be incorporated into a carcinogenesis model.

ODENames={AmtFat,AmtMuscle,AmtLiver,AmtRichly};
RatSoln=NDSolve[FinalRatODE,ODENames,{t,0,100.},MaxSteps->5000];

For simplicity, we  assume that the cellular birth rate β(t) of a (Subscript[B, 2])/0/0 model is proportional to the amount of styrene in the liver of a rat. Specifically, let β(t)=0.01+0.05×AmtLiver(t), where AmtLiver(t) is the amount of styrene in the liver of a rat predicted by the PBPK model. The survival function for the rat can now be easily computed and plotted.

RatSurvival[t_]:=B00NonhomoSurvival[t,0.1&,{{(0.01+0.05*RatSoln[[1]][[3]][[2]][#])&,0.01&,0.01&}}]

Plot[RatSurvival[s],{s,0,100},Frame->True,
PlotLabel->FontForm["A Survival Function",{"Roman",16}]];

We now recall briefly that, if S(t) denotes a survival function, the corresponding hazard  function is defined by  h(t) = -S' (t)/S(t), where S' (t) stands for  the derivative of the survival function S(t). Because the hazard  at a given time 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  at (Subscript[t, 0]). Consequently, with nonhomogeneous models, several ODE systems must be solved to get one value of the hazard  function. It is therefore natural and economical that, with time-dependent parameter models, CarcinoMod functions return interpolating function objects for hazard functions.  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. Numerical differentiation of the survival function S(t) is simple in principle but tricky in practice. Three options are available that allow the user to control how  interpolating function objects are constructed and hence allow the user to improve the numerical differentiation process. The first 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. For instance, the hazard function of the above styrene example can be drawn as follows.

styreneHaz=B00NonhomoHazard[100,0.1&,{{(0.01+0.05*RatSoln[[1]][[3]][[2]][#])&,0.01&,0.01&}},HazardInterpolationOrder->10,HazardInterpolationPoints->50,HazardBoundaryStretch->0]

Plot[styreneHaz[s], {s, 0, 100}, PlotRange -> All, Frame -> True, 
      PlotLabel -> "A Hazard Function"];

Cumulant Functions

In addition to the survival and hazard functions, cumulant functions are also useful and interesting quantities. For instance, Portier and Kopp-Schneider (1991) derived cumulant functions for certain carcinogenesis model with DNA damage and repair mechanisms. Cumulant functions are equivalent to moment functions, but are more amenable. Take a  (Subscript[B, 2])/N/0 model for example. The cumulant function (Subscript[κ, 1, 0, 0])(t) is the mean number of normal cells, while the cumulant function  (Subscript[κ, 2, 0, 0])(t) is the variance of normal cell counts. Constant parameters are assumed by CarcinoMod when computing cumulant functions. CarcinoMod returns cumulant functions as interplating function objects.

?BN0Cumulant

meanFunc = 
      BN0Cumulant[
        1000, {{0.07, 0.063, 0.0001}, {0.068, 0.062, 0.0001}, {0.07, 0.063, 
            0.0001}}, {1, 0, 0, 0}, 5]

Plot[meanFunc[s], {s, 0, 1000}, PlotRange -> All, Frame -> True, 
      PlotLabel -> "Mean Number of Normal Cells"];

We now compute the covariance function between the counts of normal cells and stage 1 cells using the same (Subscript[B, 2])/N/0 model.

v20 = 
      BN0Cumulant[
        3000, {{0.07, 0.063, 0.0001}, {0.068, 0.062, 0.0001}, {0.07, 0.063, 
            0.0001}}, {2, 0, 0, 0}, 5];

v02 = 
      BN0Cumulant[
        3000, {{0.07, 0.063, 0.0001}, {0.068, 0.062, 0.0001}, {0.07, 0.063, 
            0.0001}}, {0, 2, 0, 0}, 5];

v11 = 
      BN0Cumulant[
        3000, {{0.07, 0.063, 0.0001}, {0.068, 0.062, 0.0001}, {0.07, 0.063, 
            0.0001}}, {1, 1, 0, 0}, 5];

covar[t_] := v11[t]/Sqrt[v20[t] v02[t]]

Plot[covar[t], {t, 0, 3000}, Frame -> True, 
      PlotLabel -> "A Covariance Function"];

Sometimes the ordinary differential equations (ODEs) that are used to compute the cumulant functions are themselves of interest. CarcinoMod can generate such symbolic ODE systems for all eight types of models using the symbols specified by the user. Computational details on deriving symbolic ODEs of this type  can be found in MathSource item 0207-942 (Zheng 1995c).

?B00CumulantODE

Before proceed, we  clear all values that may be accidentally stored in the symbols  we intend to use, to ensure that parameter names, cumulant names and time variables are all treated as symbols by Mathematica.  We also require that indices be displayed as subscripts in the traditional fashion.

\(Clear[ν, β, δ, μ, κ, t]; 

    Format[κ[x__]] := Subscript[κ, x]\)\(Format[β[s_]] := Subscript[β, s]; 

    Format[δ[s_]] := Subscript[δ, s]; 

    Format[μ[s_]] := Subscript[μ, s]; 
\)

\(\(B00CumulantODE[{ν, β, δ, μ}, κ, s, 2, 
        2] // TableForm
\)\)

Note that parameters and cumulants are indexed differently between a  (Subscript[B, 2])/0/0 model and a (Subscript[B, 2])/N/0 model.

\(\(BN0CumulantODE[{β, δ, μ}, κ, s, 2, 2] // 
      TableForm
\)\)

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 functions.  In fact, all the analytic formulae derived for the survival and hazard 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". We have already encountered two of these functions: B200Survival and B200Hazard. Another function in this group is  B200PiecewiseLinearHazard which computes the hazard function for a (Subscript[B, 2])/0/0 model that allows the first mutation rate (the intensity of a Poisson process) to change with  time in a piecewise linear manner. This idea originated in a lung cancer study by Moolgavkar et al. (1989), and was also applied in a data analysis of atomic bomb survivors  by Kai et al. (1997). The formula was derived by this author in an unpublished manuscript and is useful when the tissue under consideration is still in its developing stage. Because in a (Subscript[B, 2])/0/0 model,  normal cell growth is assumed to follow some deterministic law (e.g., the Gompertz law or 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 function in a simple yet effective way. Here we assume that normal cell growth can be adequately described by a Gompertz function.

\(gompertz[a_, b_, N0_, t_] := N0\ Exp[a/b\ \((1 - Exp[\(-b\)\ t])\)]; 

    normalCellGrowth = 10\ ^\((\(-7\))\) gompertz[0.1, 0.006, 1, #]&; 

    Plot[normalCellGrowth[t], {t, 0, 1500}, Frame -> True, 
      PlotLabel -> "A Gompertz Growth Function"]; \)

?B200PiecewiseLinearHazard

It would be tedious to determine the constants (Subscript[a, i]) and (Subscript[b, i]) from a given function f and partion points (Subscript[t, i]). The function LinearizeFunction is designed to automate this process.

?LinearizeFunction

ti = Range[0, 1500, 500];

triplets = LinearizeFunction[normalCellGrowth, ti]

g1 = 
      Plot[B200PiecewiseLinearHazard[t, triplets, {0.02, 0.018, 0.001}], {t, 
          0, 1800}];

We now refine the partition of the time interval [0,1500], and compare the corresponding hazard functions.

\(ti = Range[0, 1500, 300]; 

    triplets = LinearizeFunction[normalCellGrowth, ti]; 

    g2 = Plot[
        B200PiecewiseLinearHazard[t, triplets, {0.02, 0.018, 0.001}], {t, 0, 
          1800}, PlotStyle -> RGBColor[0, 1, 0], 
        DisplayFunction -> Identity]; \)

\(ti = Range[0, 1500, 50]; 

    triplets = LinearizeFunction[normalCellGrowth, ti]; 

    g3 = Plot[
        B200PiecewiseLinearHazard[t, triplets, {0.02, 0.018, 0.001}], {t, 0, 
          1800}, PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction -> Identity]; \)

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

The last two functions in this group deal with the (Subscript[B, 2])/0/0 model whose four parameters can be piecewise constant. The algorithm for computing the survival function for this case can be found in Zheng (1995a). It was first implemented as a  Mathematica  function in MathSource item  0207-931 (Zheng 1995b).   The algorithm for computing the hazard rate function for the same case was provided by Heidenreich et al. (1997).  Survival functions are computed using B200PiecewiseConstantSurvival, and we demonstrate its use by the  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];

Hazard  functions are computed by calling  B200PiecewiseConstantHazard, and we demonstrate its use by reproducing  Figure 3 of Heidenreich et al. (1997).

<< Graphics`Graphics`

g0 = 
      LogPlot[10000 
          B200Hazard[t, 5.5,  {23.7, 23.7 *0.995, 5.5 *10^(-7)}], {
          t, 0, 120}, PlotStyle -> Dashing[{0.02}], 
        PlotRange -> {{0, 120}, {0.1, 1300}}, Frame -> True, 
        DisplayFunction -> Identity];

curve1 = {{5, 23.7, 0.995 23.7, 5.5, 5.5 10^(-7)}, {15, 
          4*23.7, 0.995 4*23.7, 5.5, 5.5 10^(-7)}, {Infinity, 23.7, 
          0.995 23.7, 5.5, 5.5 10^(-7)}};

g1 = 
      LogPlot[10000 B200PiecewiseConstantHazard[t, curve1], {t, 0, 120}, 
        PlotStyle -> RGBColor[1, 0, 0], DisplayFunction -> Identity];

curve2 = {{25, 23.7, 0.995 23.7, 5.5, 5.5 10^(-7)}, {35, 
          4*23.7, 0.995 4*23.7, 5.5, 5.5 10^(-7)}, {Infinity, 23.7, 
          0.995 23.7, 5.5, 5.5 10^(-7)}};

g2 = 
      LogPlot[10000  B200PiecewiseConstantHazard[t, curve2], {t, 0, 120}, 
        PlotStyle -> RGBColor[0, 1, 0], DisplayFunction -> Identity];

curve3 = {{65, 23.7, 0.995 23.7, 5.5, 5.5 10^(-7)}, {75, 
          4*23.7, 0.995 4*23.7, 5.5, 5.5 10^(-7)}, {Infinity, 23.7, 
          0.995 23.7, 5.5, 5.5 10^(-7)}};

g3 = 
      LogPlot[10000  B200PiecewiseConstantHazard[t, curve3], {t, 0, 120}, 
        PlotStyle -> RGBColor[0, 0, 1], DisplayFunction -> Identity];

Show[g0, g1, g2, g3, DisplayFunction -> $DisplayFunction];

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

Conolly, R.B., Reitz, R.H., Clewell, III, H.J. and Andersen, M.E. (1988). Pharmacokinetics, biochemical mechanism and mutation accumulation: a comprehensive model of chemical carcinogenesis,  Toxicology Letter, 43:189-200.

Denes, J and Krewski, D. (1996). An exact representation for the generating function for the Moolgavkar-Venzon-Knudson two-stage model of carcinogenesis with stochastic stem cell growth, Mathematical Biosciences, 131:185-204.

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.

Kai, M, Luebeck, E.G. and  Moolgavkar, S.H. (1997).  Analysis of the incidence of solid cancer among atomic bomb survivors using a two-stage model of carcinogenesis,  Radiation Research, 148:348-385.

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.

Kopp-Schneider, A. (1997). Carcinogenesis models for risk assessment, Statistical Methods in Medical Research, 6:317-340. 

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.

Moolgavkar, S.H., Dewanji, A. and Luebeck, G. (1989). Cigarette smoking and lung cancer: Reanalysis of the British doctor's data, Journal of the National Cancer Institute, 81:415-420.

Moolgavkar, S.H. and Luebeck, E.G. (1990). Two-event model for carcinogenesis: biological, mathematical and statistical considerations, Risk Analysis, 10:323-341.

Portier, C.J. and Kopp-Schneider, A. (1991). A multistage model of carcinogenesis incorporating DNA damage and repair,  Risk Analysis, 11:535-543.

Portier, C.J., Kopp-Schneider, A. and  Sherman, C.D. (1997).  Calculating tumor incidence rates in stochastic models of carcinogenesis, Mathematical Biosciences, 139:129-146.

Quinn, D.W. (1989). Calculating the hazard function and probability of tumor for cancer risk assessment when the parameters are time-dependent,  Risk Analysis, 9:407-413.
 
Ramsey, J.R. and Andersen, M.E. (1984). A physiological model for the inhalation pharmacokinetics of inhaled styrene in rats and humans,  Toxicology and Applied Pharmacology, 73:159-175.

Tan, W.Y. (1991). Stochastic Models of Carcinogenesis, Marcel Dekker, New York.

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. (1995c). MathSource item 0207-942, 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 -- Simulation and Computation, 27:137-145.