(* This package provides functions for numerical investigation of the MVK carcinogenesis models that assume Erlang distributed life lengths for initiated cells. The main reference is the paper by the same author entitled "On the MVK Stochastic Carcinogenesis Model With Erlang Distributed Cell Life Lengths" (Risk Analysis, 15:495-520, 1995). The author is indebted to David Withoff and Tom Issaevitch of Wolfram Research for providing valuable programming suggestions while he was visiting Wolfram Research during the final stage of the work. Draft version August 17, 1994, Jefferson, Arkansas Revision August 2, 1995, Champaign, Illinois *) (* Qi Zheng Division of Biometry and Risk Assessment National Center for Toxicological Research Jefferson, Arkansas 72079-9502 *) BeginPackage["CarcinoErlang`"] SurvivalMVK::usage="SurvivalMVK[b,d,nu,mu,t] computes the survival probability at time t for the classic MVK model. The parameters b,d, nu, and mu stand for the birth, death, first mutation and second mutation rates, respectively." SurvivalErlang::usage="SurvivalErlang[b,d,nu,mu,k,t0] calculates the survival probability at time t0 assuming an Erlang (order k>=2) distribution for the initiated cell life lengths. The birth, death, first and second mutation rates are b, d, nu and mu, respectively. These parameters are adjusted such that the mean cell life length is the same as in the original MVK model (k=1) with the same birth, death and second mutation rates." MeanMVK::usage="MeanMVK[b,d,nu,mu,t0] computes the mean value functions of the classic MVK model at time t0. Returned value is of the form {m1,m2}, where m1 is the mean number of initiated cells and m2 is the mean number of malignant cells. The birth, death, first and second mutation rates are b,d,nu and mu, respectively." MeanErlang::usage="MeanErlang[b,d,nu,mu,k,t0] returns a list {m1,m2}, where m1 is the mean number of the initiated cells and m2 is the mean number of the malignant cells at time t0. Life length distribution for the initiated cell is assumed to be Erlang (order k>=2). The birth, death, first and second mutation rates are b,d,nu and mu, respectively. These parameters are adjusted such that the mean cell life length is the same as in the original MVK model (k=1) with the same birth, death and second mutation rates." PlotMeanErlang::usage="PlotMeanErlang[b,d,nu,mu,k,maxT,InitMalig:=1] plots the mean value functions of the MVK model when the life time distribution of the initiated cells is assumed Erlang (order k>=2). Setting the parameter InitMalig=1 plots the mean value function for the initiated cells over the time interval (0,maxT), otherwise it plots the mean value function for the malignant cells. The birth, death, first and second mutation rates are b,d,nu and mu, respectively. These parameters are adjusted such that the mean cell life length is the same as in the original MVK model (k=1) with the same birth, death and second mutation rates." SurvivalPiecewiseConstant::usage="SurvivalPiecewiseConstant[regimen,t] computes the survival probability at time t for the MVK model whose parameters are piecewise constants. The parameter regimen is a list of lists containing the model parameters on each time interval along with the end points of the time interval. For, example, {{6000,0.0019,0.0015,0.03,3 10^(-7)}, {Infinity,0.0017,0.0014,0.02,2 10^(-7)}} partitions the time interval into (0,6000) and (6000, Infinity). The other four parameters in each sublist are birth, death, first and second mutation rates." Begin["Private`"] (* computes the survival function of the classic MVK model using the exact formula *) SurvivalMVK[b_,d_,nu_,mu_,t_]:=Module[{c1,c2,r1,r2,num,den}, c1=b+d+mu; c2=Sqrt[c1^2-4 b d]; r1=(c1-c2)/(2 b); r2=(c1+c2)/(2 b); num=(r2-r1) Exp[b (r2-1) t]; den=1-r1-(1-r2) Exp[b (r2-r1) t]; (num/den)^(nu/b) ]/;(b>=0)&&(d>=0)&&(nu>=0)&&(mu>=0)&&(t>=0) SurvivalErlang[b_,d_,nu_,mu_,k_Integer,t0_]:=Module[ {c,eqn1,eqn2,eqns,init,soln,integ,z,t}, c=b+d+mu; eqn1=Table[z[i]'[t]==k c (z[i][t]-z[i+1][t]),{i,1,k-1}]; eqn2={z[k]'[t]==-k b z[1][t]^2+k c z[k][t]-k d}; init=Table[z[i][t0]==1,{i,1,k}]; eqns=Join[eqn1,eqn2,init]; soln=NDSolve[eqns,Table[z[i][t],{i,1,k}],{t,0,t0},MaxSteps->3000]; integ=NIntegrate[Evaluate[z[1][t] /. soln[[1]] ],{t,0,t0}]; Exp[-nu t0] Exp[nu integ] ]/; (b>=0)&&(d>=0)&&(nu>=0)&&(mu>=0)&&(k>=2)&&(t0>=0) (* computes mean value functions for the classic MVK model using exact formulae*) MeanMVK[b_,d_,nu_,mu_,t0_]:=Module[ {m1,m2}, If[b==d, m1=nu t0; m2=nu mu t0^2/2, m1=(Exp[(b-d) t0]-1) nu/(b-d); m2=(m1-nu t0) mu/(b-d)]; {m1,m2} ]/;(b>=0)&&(b>=0)&&(nu>=0)&&(mu>=0)&&(t0>=0) MeanErlang[b_,d_,nu_,mu_,k_Integer,t0_]:= Module[{c,eqn1,eqn2,eqn3,eqn4,init,soln,m,t}, c=b+d+mu; eqn1={m[1]'[t]==-k c m[1][t]+k (2 b+mu) m[k][t]+nu}; eqn2=Table[m[i]'[t]==-k c m[i][t]+k c m[i-1][t],{i,2,k-1}]; eqn3={m[k]'[t]==-k c m[k][t]+k c m[k-1][t]}; eqn4={m[k+1]'[t]==k mu m[k][t]}; init=Table[m[i][0]==0,{i,1,k+1}]; eqns=Join[eqn1,eqn2,eqn3,eqn4,init]; soln=NDSolve[eqns,Table[m[i][t],{i,1,k+1}],{t,0,t0},MaxSteps->3000]; {Sum[m[i][t],{i,k}],m[k+1][t]}/.soln[[1]]/.t->t0 ]/;(b>=0)&&(d>=0)&&(nu>=0)&&(mu>=0)&&(k>=2)&&(t0>=0) PlotMeanErlang[b_,d_,nu_,mu_,k_Integer,maxT_,InitMalig_:1,opt___]:= Module[{c,eqn1,eqn2,eqn3,eqn4,init,soln,m,t}, c=b+d+mu; eqn1={m[1]'[t]==-k c m[1][t]+k (2 b+mu) m[k][t]+nu}; eqn2=Table[m[i]'[t]==-k c m[i][t]+k c m[i-1][t],{i,2,k-1}]; eqn3={m[k]'[t]==-k c m[k][t]+k c m[k-1][t]}; eqn4={m[k+1]'[t]==k mu m[k][t]}; init=Table[m[i][0]==0,{i,1,k+1}]; eqns=Join[eqn1,eqn2,eqn3,eqn4,init]; soln=NDSolve[eqns,Table[m[i][t],{i,1,k+1}],{t,0,maxT},MaxSteps->3000]; If[InitMalig==1, Plot[Evaluate[Sum[m[i][t],{i,1,k}]/.soln],{t,0,maxT},opt], Plot[Evaluate[m[k+1][t]/.soln],{t,0,maxT},opt]] ]/;(b>=0)&&(d>=0)&&(nu>=0)&&(mu>=0)&&(k>=2)&&(maxT>0) (*The zi function defined in Eq(13) with ti-t replaced by dt*) zi[ci_,{b_,d_,nu_,mu_,dt_}]:=Module[ {c1,c2,r1,r2}, c1=b+d+mu; c2=Sqrt[c1^2-4 b d]; r1=(c1-c2)/(2 b); r2=(c1+c2)/(2 b); With[{w=(r2-ci) Exp[b (r2-r1) dt]},(r2 (ci-r1)+r1 w)/(ci-r1+w)] ] (*computes a single term in Eq(18)*) S[dt_,{b_,d_,nu_,mu_},ci_]:=Module[ {c1,c2,r1,r2}, c1=b+d+mu; c2=Sqrt[c1^2-4 b d]; r1=(c1-c2)/(2 b); r2=(c1+c2)/(2 b); ((r2-r1) Exp[b (r2-1) dt]/(ci-r1+(r2-ci) Exp[b (r2-r1) dt]))^(nu/b) ] (* computes the survival probability according to Eq (18), regimen is of the form {{t1,b1,d1,mu1,nu1},...,{tm,bm,dm,mum,num}}*) PiecewiseConstant[regimen_List]:=Module[ {dtList,cList}, dtList=Rest[#-RotateRight[#]]&[First[#]&/@ regimen]; cList=FoldList[zi,1, Reverse[MapThread[Append,{Rest[Rest[#]&/@regimen],dtList}]]]; Times@@MapThread[S,{Prepend[dtList,First[First[regimen]]], Rest/@regimen,Reverse[cList]}] ] SurvivalPiecewiseConstant[regimen_List,t_]:=Module[ {nPieces,newRegimen}, nPieces=Count[First/@regimen,x_/;(x=0)||(x==Infinity)]&,#]&/@ regimen]&&((And@@((Length[#]==5)&/@ #))&[regimen]) End[] EndPackage[]