Carcinogenesis Models With Erlang Distributed Cell Life Lengths

Qi Zheng
National Center for Toxicological Research 
U.S. Food and Drug Administration

Introduction

Carcinogenesis models are an important tool in quantitative cancer research. The model that will be the subject of this notebook is based on  the classic MVK two-stage carcinogenesis model but assumes Erlang distributed cell life lengths for the initiated cells.  When the order of the Erlang distribution is one, the model reduces to the original MVK  model;  the original MVK model is hence included  for the sake of completeness. Details for the Erlang  model are contained in a paper by this author, "On the MVK Stochastic Carcinogenesis Model with Erlang Distributed Cell Life Lengths" (Risk Analysis, 15:495-502, 1995) .  The package CarcinoErlang.m consisting of 6 functions is designed to facilitate the numerical investigation reported in the  paper. This notebook  provides example uses of the package. Readers interested in the modeling aspects should consult the aforementioned paper and references therein.

<<~/mathsource/CarcinoErlang.m

Survival Functions

The package provides two functions for computing the survival functions. For the classic MVK model, the closed form analytic expression is used (see, e.g., Risk Analysis, 14:1079-1080 or 1081-1084 , 1994). For the general case, the approach is based on numerical solutions of ordinary differential equations.

?SurvivalMVK

SurvivalMVK[0.0019,0.0015,0.03,3 10^-7,9000]

?SurvivalErlang

SurvivalErlang[0.0019,0.0015,0.03,3 10^-7,5,9000]

The following example compares the survival functions when the cell life lengths of the initiated cells follow Erlang distributions of orders 1, 2 and 5. Note that  the first case  employs the function SurvivalMVK[ ].

f1=SurvivalMVK[0.0019,0.0015,0.03,0.0000003,#]&;
f2=SurvivalErlang[0.0019,0.0015,0.03,0.0000003,2,#]&;
f5=SurvivalErlang[0.0019,0.0015,0.03,0.0000003,5,#]&;
<<Graphics`Legend`

Plot[{f1[t],f2[t],f5[t]},{t,0,12000},
PlotStyle->{Dashing[{0.01}],Dashing[{0.02}],
    Dashing[{0.04}]},Frame->True,
    FrameLabel->{"Time","Survival Probability"},
    PlotLabel->FontForm["Comparing Survival Curves",
           {"Palatino-Bold",18}],
PlotLegend->{"MVK","Erlang(2)", "Erlang(5)" },
LegendPosition->{0.8,0.2},LegendSize->{0.5,0.5}]

Mean Value Functions

The package also  computes the mean value functions for both the initiated cells and the malignant cells.  Exact analytic formulae are employed for the original MVK model,  and numerical approaches are used for  other  cases. For efficiency reasons,  a separate PlotMeanErlang[ ] function  plots the mean value functions for the Erlang (order k>=2)  cases.

?MeanMVK

MeanMVK[0.0019,0.0015,0.03,3 10^-7,7500]

?MeanErlang

MeanErlang[0.0019,0.0015,0.03,3 10^-7,5,7500]

To draw the mean value function for the initiated cells for the original MVK model, it may be convenient to define a temporary  function .

meanInit=(MeanMVK[0.0019,0.0015,0.03,3 10^-7,#][[1]])&;
Plot[meanInit[t],{t,0,12000},Frame->True,
FrameLabel->{"Time","Mean Number"},
PlotLabel->FontForm["Mean Number of Initiated Cells",
{"Palatino-Bold",18}]]

The following  compares  the mean value functions for the malignant cells when the life length distributions for the initiated cells follow Erlang distributions of orders 1, 2 and 5.

?PlotMeanErlang

malig1=(MeanMVK[0.0019,0.0015,0.03,3 10^-7,#][[2]])&;
p1=Plot[malig1[t],{t,0,12000},Frame->True,PlotRange->All,
DisplayFunction->Identity,
FrameLabel->{"Time","Mean Number"},
PlotLabel->FontForm["Mean Number of Malignant Cells",
{"Palatino-Bold",18}]];
p2=PlotMeanErlang[0.0019,0.0015,0.03,3 10^-7,2,12000,2,
PlotStyle->{Dashing[{0.01}]},
DisplayFunction->Identity];
p5=PlotMeanErlang[0.0019,0.0015,0.03,3 10^-7,5,12000,2,
PlotStyle->{Dashing[{0.03}]},
DisplayFunction->Identity];
Show[p1,p2,p5,DisplayFunction->$DisplayFunction]

The MVK Model With Piecewise Constant Parameters

The package includes a function  calculating the survival probabilities for the MVK model when the model parameters are  piecewise constants.

?SurvivalPiecewiseConstant

The following shows an imaginary situation with time  divided into five subintervals.  Note the end point of the last time interval is always the built-in symbol Infinity.

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}};

What does the survival function look like for this complicated case?

surv=SurvivalPiecewiseConstant[fivePieces,#]&;
Plot[surv[t],{t,0,20000},Frame->True,
FrameLabel->{"Time","Survival Probability"},
PlotLabel->FontForm["A Survival Curve",
{"Palatino-Bold",18}],GridLines->Automatic];

A  graph similar to Figure 5 in the main reference is reproduced here. This type of  cases may be interpreted as  different dose patterns in an animal experiment.

s=SurvivalPiecewiseConstant[
   {{#1,0.0019,0.0015,0.03,3 10^(-7)},
    {Infinity,0.0015,0.0015,0.03,3 10^(-7)}},#2]&;
Plot[{s[1000,t],s[3000,t],s[6000,t],s[16000,t]},{t,0,25000},
PlotStyle->{Dashing[{0.01}],Dashing[{0.02}],
   Dashing[{0.01,0.03,0.01}],Dashing[{0.04}]},
LegendLabel->"Changing\n Point",
PlotLegend->{"1000","3000","6000","16000"},
LegendPosition->{0.7,0.1},LegendSize->{0.5,0.6},
Frame->True,FrameLabel->{"Time","Survival Probability"},
PlotLabel->FontForm["Four Survival Curves",
{"Palatino-Bold",18}]];

Acknowledgments

A  majority of this work was performed while the author was visiting Wolfram Research, Inc.. I am grateful to Wolfram Research for providing a hospital atmosphere and partial financial support for completing  this research.  The author is indebted to David Withoff  and Tom  Issaevitch of Wolfram Research for providing valuable programming suggestions. This work was also  supported by an appointment to the Postgraduate Research Program at the National Center for Toxicological Research administered by the Oak Ridge Institute for Science and Education through an interagency agreement between the U.S. Department of Energy and the Food and Drug Administration.