Cells mutation: a discrete-time Markov chain model​
​By Alessandro Esposito​
​The Medical Research Council Cancer Unit,​
​ University of Cambridge, Hills Road, Cambridge CB2 0XZ, UK

Abstract

Most tumours exhibit significant heterogeneity and are best described as communities of cellular populations competing for resources. Growing experimental evidence also suggests that cooperation between cancer clones is important as well for the maintenance of tumour heterogeneity and tumour progression. However, a role for cell communication during the earliest steps in oncogenesis is not well characterized despite its vital importance in normal tissue and clinically manifest tumours. Here, we present a simple analytical model and stochastic lattice-based simulations to study how the interaction between the mutational process and cell-to-cell communication in three-dimensional tissue architecture might contribute to shape early oncogenesis. We show that non-cell-autonomous mechanisms of carcinogenesis could support and accelerate pre-cancerous clonal expansion through the cooperation of different, non- or partially transformed mutants. We predict the existence of a ‘cell-autonomous time horizon’, a time before which cooperation between cell-to-cell communication and DNA mutations might be one of the most fundamental forces shaping the early stages of oncogenesis. The understanding of this process could shed new light on the mechanisms leading to clinically manifest cancers.

Two-hits model

We define a discrete-time Markov chain model based on four states: W, X, Y, and XY. W represents the wild-type state of a cell, X and Y are cell states characterized by two independent mutations (X and Y), and XY the cell that have accrued both mutations. To compare the Markov chain model with other calculations provided in our work, we assume that cells are subject to a low (p<<1 mutation/day) and constant mutational rate. The transition matrix T can be then described as:
In[]:=
T=
1-xp-yp-xy
2
p
xp
yp
xy
2
p
0
1-yp
0
yp
0
0
1-xp
xp
0
0
0
1
;
In[]:=
The initial state vector (S0) will represent the system where every cell is initially wild-type.
In[]:=
S0=(
1
0
0
0
);
We can then define the discrete-time Markov process with  and the probability density function for the observation of the sate 4 (i.e., XY) after k days as:
In[]:=
=DiscreteMarkovProcess[S0[[1]],T];​​pxy[k_]=PDF[[k],4]
Out[]=
1-
k
(1-px)
1+px
-
k
(1-py)
1+py
+
(1-
2
p
xy)
k
(1-px-py-
2
p
xy)
(1+px)(1+py)
As p is small by definition, we can approximate
p
xy
(k) with the first non-zero term of a Taylor expansion.
In[]:=
pxy[k_]=Normal[Series[PDF[[k],4],{p,0,2}]]
Out[]=
2
k
2
p
xy
Within a set of N independent cells, the probability that none of them is mutant will be
N
(1-
p
xy
(k))
and the probability to observe a mutant cell at time k is thus: (note that we used the symbol of capital “Nu” to represent the letter “N” in the code to prevent contradiction with the function N[ ])
In[]:=
Pxy[k_]=1-
Ν
(1-pxy[k])
Out[]=
1-
Ν
(1-
2
k
2
p
xy)
For
2
k
2
p
xy
small (i.e. for k not to large as expected during the early stages of the modelled process of carcinogenesis)
1-
p
xy
(k)
~ 1-
2
k
2
p
xy
~
2
k
2
p
xy
e
permitting us to describe
P
xy
in the rather convenient form of
In[]:=
Pxy[k_]=1-Exp[-Νxy
2
k
2
p
];
P
xy
(k)
is the cumulative distribution function related to observe a double XY mutant cells in a population of N cells accruing X and Y mutations with a probability of xp and yp, respectively. The PDF can be thus evaluated by differentiating
P
xy
(k)
.
In[]:=
PDFxy[k_]=D[Pxy[k],{k,1}]
Out[]=
2
-
2
k
2
p
xyΝ

k
2
p
xyΝ
Finally this PDF can be used to evaluate the expectation of the mean first passage time, i.e. the time at which on average we would observe the first XY mutant cells in a population of N cells.
In[]:=
km=Integrate[kPDFxy[k],{k,0,∞},Assumptions{Ν>1,p>0,p<1,Ω>0,x>0,y>0}]
Out[]=
π
2p
xyΝ
k
m
for two independent mutations occurring in a single cells (x=y=1) will be thus
In[]:=
ta=km/.{x1,y1}
Out[]=
π
2p
Ν
k
m
for two independent mutations occurring in two distinct cells within a neighbourhood defined by the parameter Ω (x=1 and y=Ω) will be
In[]:=
tΩ=km/.{x1,yΩ}
Out[]=
π
2p
ΝΩ
The scaling factor
t
Ω
/
t
a
is thus
In[]:=
Simplify
tΩ
ta
,AssumptionsΝ>0
Out[]=
1
Ω
In the published manuscript, this model is used to compare analytical results with Monte Carlo simulations (Fig. 3d in https://doi.org/10.1098/rsos.201532). Here, we show how the distribution of the waiting time for CD clones at different values of Ω. For two mutations to occur within the same cells (equivalent to Ω=1), the distributions of the waiting times is peaked at
t
a
and is very broad. For larger values of the oncogenic field effect, i.e. when mutation can occur in different cells within a given neighbourhood but mediating a pre-cancerous phenotype by cell-to-cell communication, the distribution is rescaled.
​
Ωlist={1,2,4,8,16,24};​​​​Show​​TableLogLinearPlotEvaluate[(12/365)*PDFxy[k]/.{Ν10^6,p10^-6,x1,yΩlist[[i]]}],{k,0,6*365},​​PlotRangeAll,PlotLabelsCalloutStringJoin["Ω=",ToString[Ωlist[[i]]]],√π2p√(Νy)/.{Ν10^6,p10^-6,x1,yΩlist[[i]]},Above,Appearance"CurvedLeader",FillingAxis,FillingStyleDirective[Opacity[.3],RGBColor[1-1/i,0,1/i]],AxesLabel{"Waiting time (simulation days)","Probability (​
-1
day
​)"},PlotStyleDirective[Opacity[1],Black,Dashing[{}],Thin],{i,1,6}​​
Out[]=

Three-hits model

In this section we repeat the same steps shown for the two-hits model but for three independent mutations, x, y, and z.
In[]:=
T=
1-xp-yp-zp-xy
2
p
-xz
2
p
-yz
2
p
-xyz
3
p
xp
yp
zp
xy
2
p
xz
2
p
yz
2
p
xyz
3
p
0
1-yp-zp-yz
2
p
0
0
yp
zp
0
yz
2
p
0
0
1-xp-zp-xz
2
p
0
xp
0
zp
xz
2
p
0
0
0
1-xp-yp-xy
2
p
0
xp
yp
xy
2
p
0
0
0
0
1-zp
0
0
zp
0
0
0
0
0
1-yp
0
yp
0
0
0
0
0
0
1-xp
xp
0
0
0
0
0
0
0
1
;​​S0=(
1
0
0
0
0
0
0
0
);​​=DiscreteMarkovProcess[S0[[1]],T];​​pxyz[k_]=PDF[[k],8]
Out[]=
k
(1-px-py-
2
p
xy-pz-
2
p
xz-
2
p
yz-
3
p
xyz)
(-
2
x
y-x
2
y
-p
2
x
2
y
+
2
p
3
x
2
y
+
⋯123⋯
)
2
(1+px)
2
(1+py)
2
(1+pz)
(y+z+pyz)(
2
x
+xy+p
2
x
y+xz+p
2
x
z+yz+2pxyz+
2
p
2
x
yz)
+
⋯1⋯
⋯1⋯
+
⋯4⋯
+
⋯1⋯
⋯1⋯
+
(-1+
2
p
xy)
⋯4⋯
(1+py)
⋯3⋯
(
2
x
y+
⋯107⋯
+
9
p
4
x
4
y
4
z
)
large output
show less
show more
show all
set size limit...
Which can be approximated with the first factor of a Taylor expansion that in the case of 3 independent mutations will be the third-order term.
pxyz[k_]=Normal[Series[PDF[[k],8],{p,0,3}]];​​Pxyz[k_]=1-Exp[-Νpxyz[k]];​​PDFxyz[k_]=D[Pxyz[k],{k,1}];​​km=Integrate[kPDFxyz[k],{k,0,∞},Assumptions{Ν>1,p>0,p<1,Ω>0,x>0,y>0,z>0}];​​ta=km/.{x1,y1,z1};​​tΩ=km/.{x1,y1,zΩ};​​Simplify
tΩ
ta
,AssumptionsΝ>0
Out[]=
1
1/3
Ω
Examples of triple clones with two mutations acting through paracrine effects