A Carcinogenesis Model  
at the Nucleotide Level

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

Introduction

Rapid and accelerating advances in biomedical cancer research have led a dramatic increase in the amount of attention given to mathematical modeling of carcinogenesis. The advent of comprehensive symbolic computation systems like Mathematica has facilitated  research  in the area. A cancer model that describes mutational events at the nucleotide level has been proposed recently by Q. Zheng, D. Gaylor and W. Lutz (American Statistical Association 1995 Proceedings of the Biometrics Section, pp. 349-354). The computational aspects of this model (both symbolic and numerical)  which were appreciably more complex than those of the traditional two-stage model   were investigated using Mathematica and reported by this author at the Tenth International Conference on Mathematical and Computer Modelling and Scientific Computing in Boston on July 6, 1995 (see "Combining Symbolic and Numerical computing in Carcinogenesis Modeling" in the Proceedings of the Tenth International Conference) .  The present  notebook is based on the results presented at the Conference and focuses on computational aspects. For biological and modeling aspects, the interested reader should consult  the two aforementioned references.

SetOptions["stdout",PageWidth->70];

Cumulants and the Cumulant Generating Function

How to Derive  Differential Equations

The model is essentially a four-compartment  stochastic process. The joint cumulants of the process provide valuable statistical information.  To solve for these cumulants,  first derive a system of ordinary differential equations (ODE) that the cumulants of up to a given order must satisfy. Such ODEs may be derived from a partial differential equation (PDE) that the cumulant generating function (cgf) satisfies. The  K[ ] function defined below  produces symbolically  the cgf for the model truncated at a specified order ord. Each term k[i1,i2,i3,i4][t]  represents a cumulant function of order i1+i2+i3+i4.

K[ord_]:=Select[ 
Sum[k[i1,i2,i3,i4][t] x1^i1 x2^i2 x3^i3 x4^i4 /(i1! i2! i3! i4!), 
   {i1,0,ord},{i2,0,ord},{i3,0,ord},{i4,0,ord}]/.k[0,0,0,0][t]->0, 
   ((Exponent[#,x1]+Exponent[#,x2]+Exponent[#,x3]+ 
   Exponent[#,x4])<=ord)&];

The right hand side of the aforementioned PDE  should also be expanded to the same  given order.

RHS[ord_]:=Select[ Expand[(b Sum[x1^i/i!,{i,1,ord}]+ 
      a Sum[(-x1+x2)^i/i!,{i,1,ord}]+ 
     d1 Sum[(-x1)^i/i!,{i,1,ord}]) D[K[ord],x1]+ 
     (p1 Sum[(x1-x2)^i/i!,{i,1,ord}]+p2 Sum[(x1-x2+x3)^i/i!,{i,1,ord}]+ 
     d2 Sum[(-x2)^i/i!,{i,1,ord}]) D[K[ord],x2]+ 
    (p3 Sum[(-x3+x4)^i/i!,{i,1,ord}]+ 
    u Sum[x4^i/i!,{i,1,ord}]+ 
    d3 Sum[(-x3)^i/i!,{i,1,ord}]) D[K[ord],x3]+ 
    v Sum[x1^i/i!,{i,1,ord}] ] , 
    ((Exponent[#,x1]+Exponent[#,x2]+Exponent[#,x3]+Exponent[#,x4])<=ord)&]

By equating coefficients of the left hand side of the PDE (which is simply D[K,t])  with those of the right hand side, a system of ODEs results.

odeForCumulants[ord_]:=List@@LogicalExpand[ 
     Series[D[K[ord],t],{x1,0,ord},{x2,0,ord},{x3,0,ord}, 
     {x4,0,ord}]==Series[RHS[ord], 
     {x1,0,ord},{x2,0,ord},{x3,0,ord},{x4,0,ord}]];

The Mean and Variance-Covariance Functions

The first and second order cumulants are the best known and the most widely used. The first order cumulants are referred to as means or expectations,  while the second order cumulants are called variances (or covariances). One can easily obtain an ODE system for these cumulants.

odeOrderTwo=odeForCumulants[2]; 
odeOrderTwo[[5]]

As one may readily see from the above output,  the resulting  ODEs are  not  yet  in a form pleasing to the human eye. Further symbolic processing is  needed to render the equations to a  "canonical" form. First, fractions should be eliminated. Because the equations are  in the form  lhs==0, one only needs to process the left hand sides of  these equations.

odeOrderTwo=(Numerator[Together[#]])&/@(First/@odeOrderTwo); 
odeOrderTwo[[5]]

Secondly, collect coefficients for each cumulant .

rule1=(a_ k[i__][t_]+ b_ k[i__][t_]):>(a+b) k[i][t]; 
odeOrderTwo=odeOrderTwo//.rule1; 
odeOrderTwo[[5]]

Thirdly, after locating the derivative term in each expression, translate these expressions back to  ODEs.

posOfDerivative=Flatten[(Position[#,_. k[__]'[_]])&/@odeOrderTwo]; 
odeOrderTwo=MapThread[(#[[#2]]==-Drop[#,{#2}])&, 
                                          {odeOrderTwo,posOfDerivative}]; 
odeOrderTwo[[5]]

Further simplification and sorting transform the ODEs into textbook format. At this point, one may  rename these cumulants with more common notation (e.g., m stands for mean, etc.).

rule2={k[1,0,0,0]->m1,k[0,1,0,0]->m2, 
k[0,0,1,0]->m3,k[0,0,0,1]->m4,k[2,0,0,0]->v11, 
k[1,1,0,0]->v12,k[1,0,1,0]->v13,k[1,0,0,1]->v14, 
k[0,2,0,0]->v22,k[0,1,1,0]->v23,k[0,1,0,1]->v24, 
k[0,0,2,0]->v33,k[0,0,1,1]->v34,k[0,0,0,2]->v44}; 
 
TableForm[Sort[Simplify[odeOrderTwo//.rule2]]]

How to Obtain Numerical Solutions

There are 34 cumulants of order up to 3. To numerically solve and plot them require  additional symbolic manipulation.  For the present purpose, one  needs these ODEs in  their "raw" form.

odeOrderThree=odeForCumulants[3]//.rule1;

Now specify numerical values for the cell kinetic parameters.

odeOrderThree=odeOrderThree//.{v->0.03,b->0.0019,d1->0.0015, 
    d2->0.0016,d3->0.0017,a->0.0005, 
    p1->0.0005,p2->0.0005,p3->0.0005,u->0.0000003};

In addition to creating a list of equation names, determine and attach  the initial conditions to the original ODE list.

odeNames=Table[If[i1+i2+i3+i4<=3,k[i1,i2,i3,i4],{}], 
                       {i1,0,3},{i2,0,3},{i3,0,3},{i4,0,3}]//Flatten//Rest; 
initial=Table[If[i1+i2+i3+i4<=3,k[i1,i2,i3,i4][0]==0,{}], 
                      {i1,0,3},{i2,0,3},{i3,0,3},{i4,0,3}]//Flatten//Rest; 
odeWithInit=Join[odeOrderThree,initial];

Note that the NDSolve function in Mathematica Version 2.2 does not always behave correctly when the functions have indexed names, like k[2,0,0,1].  When computing numerical values, NDSolve sometimes also computes numerical values for the indices, after which it is no longer able to recognize the functions.  This limitation of NDSolve has been corrected for Version 3 of Mathematica.  However, there are several simple ways to deal with this problem in Version 2.2. One of the easiest ways is to use the NProtectedAll attribute to protect the function indices from numerical computation.  The NProtectedAll attribute was added to Mathematica after the manual was printed, and so is not included in the manual.  The name  NProtectedAll has been changed to NHoldAll for Version 3.

SetAttributes[k, NProtectedAll]; 
odeSoln=NDSolve[odeWithInit,odeNames,{t,0,5000}];

To have a feeling of the shapes of these cumulant functions,  one can plot  a few of them .

pic=(Plot[Evaluate[#[t]/.odeSoln],{t,0,5000}, 
                Frame->True,PlotRange->All,PlotLabel->ToString[#], 
                DisplayFunction->Identity ] )&/@odeNames[[Range[27,34]]]; 
Show[GraphicsArray[Partition[pic,2]]]

Survival and Hazard Functions

The survival and hazard functions play a key role in carcinogenesis modeling. They help answer a variety of questions arising from quantitative cancer research. This section  briefly demonstrates computing and plotting of  these two important functions.  
 
The survival  function is computed using the characteristic method from  standard partial differential equation  theory.

Survival[a_,b_,d1_,d2_,d3_,p1_,p2_,p3_,u_,v_,t0_]:=Module[ 
  {ode,init,odeSoln,z1,z2,z3,t}, 
  ode={z1'[t]==-b[t] z1[t]^2+(a[t]+b[t]+d1[t]) z1[t]-a[t] z2[t]-d1[t], 
       z2'[t]==-p2[t] z1[t] z3[t]-p1[t] z1[t]+(p1[t]+p2[t]+d2[t]) z2[t]-d2[t],  
       z3'[t]==(u[t]+p3[t]+d3[t]) z3[t]-d3[t]}; 
  init={z1[t0]==z2[t0]==z3[t0]==1}; 
  odeSoln=NDSolve[Join[ode,init],{z1,z2,z3},{t,0,t0},MaxSteps->3000]; 
  NIntegrate[Evaluate[v[t] (z1[t]-1) /.odeSoln[[1]] ],{t,0,t0}, 
                  MaxRecursion->10]//Exp];

Note that  this approach allows the cell kinetic parameters to  be time-dependent. The following example computes the survival probabilities using  a linear function for  one of the  model parameters.

survivalList=Table[ 
   {i,Survival[0.0005&,0.0019&,0.0015&,0.0016&,0.0017&,0.0005&, 
    0.0005&,0.0005&,0.0000003&,0.03 (1+0.0000015 #) &,i]}, {i,0,6000,100}];

ListPlot[survivalList,PlotJoined->True,Frame->True, 
PlotLabel->FontForm["A Survival Curve",{"Palatino-Bold",18}], 
 FrameLabel->{"Time ", "Survival Probability"}]

The hazard function is computed by first creating an interpolation function object so that derivatives of the survival function may be computed.

interpObj=Interpolation[survivalList,InterpolationOrder->5]; 
dsurvival=Table[Evaluate[interpObj'[t]],{t,0,6000,100}];

From the definition of the hazard function, one can easily compute and plot the hazard function.

hazardList=-dsurvival/(#[[2]]&/@survivalList); 
 
ListPlot[hazardList,PlotJoined->True,Frame->True, 
 FrameLabel->{"Time  * 100","Hazard Rate"}, 
 PlotLabel->FontForm["A Hazard Curve",{"Palatino-Bold",18}]]

Acknowledgments

Part of this work was performed while the author was visiting Wolfram Research, Inc.. During this visit, Tom Issaevitch and David Withoff of Wolfram Research made several valuable and specific coding suggestions, which led to more elegant and efficient code  in the present  notebook.   I am  indebted to S.R. Sipcic at Boston University College of Engineering for kindly inviting me to present part of this work to the Tenth International Conference on Mathematical and Computer Modeling and Scientific Computing. Partial  support provided by Wolfram Research during the visit is gratefully acknowledged. 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.