A one-dimensional parameter-free model for carcinogenesis in gene expression space

A small portion of a tissue defines a microstate in gene expression space. Mutations, epigenetic events or external factors cause microstate displacements which are modeled by combining small independent gene expression variations and large Levy jumps, resulting from the collective variations of a set of genes. The risk of cancer in a tissue is estimated as the microstate probability to transit from the normal to the tumor region in gene expression space. The formula coming from the contribution of large Levy jumps seems to provide a qualitatively correct description of the lifetime risk of cancer in 8 tissues, and reveals an interesting connection between the risk and the way the tissue is protected against infections.

www.nature.com/scientificreports/ Recall the interpretation of points in Fig. 1 top panel. Each point comes from a small sample, the GE data obtained from it contains the contribution of many cells and the complex signaling system regulating their interactions. One may speak of a tissue microstate. On the other hand, points come from different patients, each carrying a particular genetic load. The fact that the points are grouped in definite regions means that these regions are indeed attractors in GE space.
We want to describe the genesis of a tumor, that is the time evolution of a portion or sample of a tissue that starts in the normal region and progress towards the tumor zone. We have already defined a single coordinate describing this progression: x 1 . In order to proceed further with the model, we shall clarify why and how this progression takes place.
The coordinate x 1 describing the tissue microstate starts at a point near the origin and realizes random oscillations in the normal zone. The cause for such random displacements is discussed in the next sections. The motion is confined to the normal region for a long time because this zone is a local maximum of fitness 10,13 . We have schematically represented in Fig. 1 bottom panel the fitness distribution along the PC1 axis. The y axis of this figure is the fitness with a minus sign, thus that the normal and tumor zones are local maxima of fitness. In the figure, the fitness is estimated from the histogram of samples along PC1. We have computed in Ref. 13 the number of available microstates in each zone, showing that this number is much greater for tumors than for normal states. In other words, the volume of the basin of attraction is much greater in the tumor than in the normal region. In addition, as a consequence of breaking the restrictions imposed by homeostasis, the mitotic rate of tumor stem cells is usually greater than that of normal somatic stem cells 14 . The conclusion is that the tumor minimum should be the deepest in Fig. 1, the one with highest fitness. Our drawing for the fitness distribution is a sketch built from the available data, however we are convinced that it is a qualitatively correct representation of the actual fitness distribution.
The intermediate region, R n < x 1 <x 1 − R t , holds a low-fitness barrier 10,13 , as shown in Fig. 1 bottom panel, which prevents the spontaneous transitions from the normal to the tumor region. The relative scarcity of samples in this region evidences the existence of the barrier.
A tissue microstate realizes random displacements within the normal region. Only when the barrier is surpassed and the microstate leaves the normal basin of attraction it is driven towards the tumor attractor. The transition is seen as discontinuous 10 .
A precise description of the transition requires the detailed knowledge of the fitness landscape and the causes of the random fluctuations. However, in order to estimate the risk of cancer in a tissue we may proceed in a simpler way and compute the probability for the x 1 variable describing the microstate to transit from the normal to the tumor region. The minimal walk length is x 1 − R n − R t . This is the goal we are aimed at in the present paper.
The starting point in our model is a large set of samples or microstates located near the origin of Fig. 1. They represent small portions of the healthy tissue. We may think of colon crypts in the studied example. The mean number of crypts in a healthy individual is estimated in Ref. 15 as 1.5 × 10 7 . We shall follow the random oscillations in GE space of each of these crypts.
With regard to the time variable, it is natural to follow the renewal cycle of somatic stem cells, guaranteeing crypt homeostasis. In the studied example, the renewal rate is 73 per year 16 . Thus, we shall measure time in terms of somatic stem cell generations. t = 0 may refer to conception or to the moment at which the first colon stem cell appears. On the other hand, t 0 = log 2 N sc , where N sc is the number of stem cells in the tissue, is the moment at which the tissue is formed. In colon, N sc ≈ 2 × 10 8 , and t 0 ≈ 27 . This is our starting point.

Small random displacements in GE space
Any variation in the expression of genes is a displacement in GE space. We conceptualize two kinds of GE variations: small displacements and large rearrangements. Naively, one may relate small displacements to variations in the expression of one or a few genes, whereas large GE rearrangements are coordinated variations of the expressions of many genes.
Small variations of GE levels spontaneously occurs and may have different origins. First, somatic mutations in the human genome are known to occur at a rate of 8 per cell generation 17 . Second, there is also a rate of accumulation of epigenetic (mainly methylation and phosphorylation) events modifying the normal expression levels 18 . Both processes could be boosted by inherited mutations 19,20 or external carcinogens 21 .
We may thus write for the x 1 coordinate, characterizing the microstate of a crypt at time t = n + 1 , the following equation: where and δê i corresponds to a random variation of the expression of the i-th gene. Eq. (1) describes a Markov chain of events 22 . On the other hand, Eq. (2) shows that fluctuations in the expression levels are filtered by the v 1 vector.
In Fig. 2 we draw the 30 genes with the greatest contributions to v 1 in COAD 11 . Positive, v 1i > 0 , and negative, v 1i < 0 , amplitudes correspond respectively to over-and under-expressed (silenced) genes in the tumor progression. We have distinguished the genes CST1 and AQP8. The former is a known marker of colon cancer 23 , whereas the latter plays a significant role in colon homeostasis 24 and should be silenced in tumors.
The maximum value of |v 1i | defines a scale, D, for the fluctuations of x 1 . In COAD, it coincides with the modulus of the v 1i related to the AQP8 gene. In order to get a simple estimate for the cancer risk, we may adopt the www.nature.com/scientificreports/ following model for the fluctuations: δx 1 = D r , where r is a uniformly distributed random number in (-1,1). This model may result from an independent variation hypothesis, i.e. random amplitudes and signs in the individual gene variations δê i , so that most of them cancel out. In this way, Eq. (1) for the small displacements in GE space describes a 1D Brownian or Poisson process 25 . We may use the well known fact that in a Brownian process, the final amplitudes at a given time are normally distributed, i.e. the probability density is given by: where a = 2/(D 2 t) . We shall evaluate the probability for a trajectory starting in the normal zone to reach the tumor zone. Above, we pointed out that the minimal walk length is R =x 1 − R n − R t . Thus, an estimate for the risk may be obtained from: where Erfc(z) is the complementary error function. The argument of this function is z = √ aR 2 = √ 2/t R/D , in principle a large number. Then, we may use the asymptotic behavior Erfc(z) ≈ exp(−z 2 )/( √ πz) for large z. The risk of cancer in COAD is obtained by multiplying the escape probability for a single crypt by the number of crypts, or by the number of stem cells, which is proportional to it: or This expression is general enough to be applied to other tissues, besides colon. The constant in Eq. (6) may account for other effects as, for example, the role of the immune system. Microregions escaping the normal region and forming a prototumor could be the subject of an attack by the immune system in the very early stages 26 . By definition, the constant is less than zero because the overall constant in Eq. (5) is less than one.
In Table 2 we compile a set of parameters for a group of tumors. The geometry of the normal and tumor regions, i.e. the parameters x 1 , R n and R t come from Ref. 10 . The D value is estimated as the maximum of |v 1i | 11 . On the other hand, the number of tissue stem cells, N sc , the stem cell turnover rate, m sc , and the lifetime risk of cancer (when available) are borrowed from Refs. 27,28 . The reported values of risk represent averages over 380 cancer registries from different cities and countries around the world 28 .
We may test Eq. (6) for the risk of cancer in a tissue resulting from small random variations of GE levels by using the data included in Table 2. A plot of the l.h.s. vs the r.h.s. of Eq. (6) should lead to a straight line with a slope near one and a constant less than zero. Notice that the life expectancy in Ref. 27 is assumed to be 80 years. Thus, t is obtained by multiplying the stem cell rate, m sc , by 80 years.
The results of that test are shown in Fig. 3. We get a nearly flat curve (slope = 2.1 × 10 −5 , or 1.5 × 10 −4 if we leave LUAD and THCA out of the fit), indicating that the proposed dependence of the risk on the parameters is not correct. Thus, the observed risk of cancer can not be explained by random variations of small amplitude in GE values. In the next section, we shall consider large GE rearrangements, or equivalently large jumps in GE space.
Let us stress that we use an expression like t = m sc × age in a very broad age interval. It is well known that m sc experiences a significant decrease as a result of aging 29,30 . However, also as a consequence of aging there is an

Large (Levy) jumps in GE space
Besides small random displacements, related to quasi independent variations in the GE values, there is also the possibility of large jumps in GE space. The origin of such large motions could be diverse. First, there are large scale mutations, involving DNA rearrangements and simultaneously modifying the expression of many genes. An example, known to play an important role in cancer, is that of aneuploidies 31 .
Second, large jumps in GE space could be related to coordinated variations in a group of genes. Indeed, GE values are known to be regulated by GE networks 32 . The global states of these networks define attractors 5,6 . Table 2. A set of parameters compiled for a group of tumors. The geometry of the normal and tumor regions, i.e. the parameters x 1 , R n and R t come from Ref. 10 . The minimal distance between both regions is R =x 1 − R n − R t . The D value is estimated as the maximum of |v 1i | 11 . On the other hand, the number of tissue stem cells, N sc and the stem cell turnover rate, m sc , are borrowed from Refs. 27,28 . The lifetime risk of cancer and its deviation (when available) is computed from Ref. 28 as the mean value and the standard deviation of the cumulative risk at a maximum age of 80 years. Bold marked tissues correspond to cancer types for which all the data is available.   Table 2 is used to this end. A very small slope is obtained in both, the full fit and a fit without LUAD and THCA, thus small amplitude fluctuations in gene expression space may not account for the risk of cancer in these tissues. www.nature.com/scientificreports/ Variations in genes playing a decisive role in the network, or accumulation of variations in many genes, may cause a transition from one of these global states to another one. Third, there is also the possibility of a programmed chain of GE variations leading to cancer, triggered by yet unknown causes, which is the basic hypothesis in the atavistic theory of cancer 33 .
For the large GE variations, we shall specify their rate of occurrence, µ , and the probability distribution for their amplitudes, π(�x 1 ) . It is very plausible to assume that π is of Pareto 34 or Levy 35 kind, with a power-like tail. Indeed, the Pareto character of GE distribution functions was demonstrated in Ref. 36 (see also 10 ). The Levy character of the length distribution functions in mutations was shown in 37 .
Thus, our assumption is that displacements in GE space are a kind of Levy flights. Small variations allow the exploration of the fitness landscape at lower scales, whereas sporadic large jumps allow to find global maxima. Besides mutations 37 , Levy flights are known to take place in many other biological processes, for example foraging 38 .
For lage | x 1 | , the tail of π is described by a Pareto exponent ν: The probability of a large jump reaching the tumor region is thus proportional to and the risk of cancer in a tissue: where we assume ν > 1 . Below, we use ν = 2 in order to get an estimate of the risk. Let us examine Eq. (9) in more details. First, Eq. (8) assumes that R is in the tail of the distribution function. This is justified if we compare R with the scale D, that in the case of COAD take the values of R = 115.65 and D = 0.0526 respectively. Second, no more than one hit or large jump is assumed to occur in the evolution of each microstate. In other words, the probability µt is less than one and large jumps are thought to be rare. Third, we should consider the possibility of large GE variations in the development period, that is why we included t 0 = log 2 (N sc ) in the formula. This is particularly important in tissues with slow renewal rates but large number of stem cells. For example, in lung t 0 ≈ 30 , but m sc × 80 years is only 5.6. Fourth, the rate of large jumps, µ , is unknown. However, if we assume roughly the same value for all tissues, then it can be absorbed in the overall constant entering Eq. (9). The Pareto exponent is also unknown. Notice that in the GE distribution functions of COAD the exponents take values between 1.6 and 2.0 10 . The value we use for estimates, ν = 2 , is motivated by this result.
Finally, we get the following expression for the risk, which may be tested against the data in Table 2: where t = t 0 + m sc × age . The constant should be negative according to our hypothesis of µ small. The results of the test with an average life expectancy of 80 years are shown in Fig. 4. The observed behavior is consistent with a linear dependence with slope near one. Indeed, we obtain a slope equal to 0.82. The Pearson correlation coefficient is 0.85, indicating that 72% of the dispersion of points may be explained by a linear dependence. The p-value is equal to 0.04. The small error bars suggest that the main reason (7) π(�x 1 ) ∼ 1/|�x 1 | ν .   Table 2. The slope of the linear fit is near one, as expected. 72 % of the data dispersion is explained by the linear dependence (p-value = 0.04). www.nature.com/scientificreports/ for the unexplained dispersion of points could be the assumption that the rate of large jumps, µ , is roughly the same for all of the tumors. A tissue specific µ variable would account for the dispersion.
In conclusion, we get the following simple expression for the risk of cancer per stem cell in a tissue, coming from large jumps in GE space: where we included an effective rate, µ ′ . Genetic, viral or external carcinogenic factors may increase µ ′ , whereas the action of the immune system in the tissue may modify µ ′ in any direction. In the next section, we qualitatively analyze a larger set of tumors by using Eq. (11).

Qualitative analysis of the data on cancer risk in different tissues
We use the extremely simple expression for the cancer risk in a tissue, offered by Eq. (11), in order to re-examine the data presented in paper 27  This expression provides a simple explanation for the intuitive claim in Ref. 27 that the risk is related to the number of stem cell replications. The results are shown in Fig. 5 and Table 3. We should try to understand the observed values of ERS in terms of the tissue characteristics. In order to facilitate the analysis, the studied tumors are separated in groups.
Group I includes 12 tumors (10 tissues), located in a band delimited by red dashed lines in Fig. 5, and coefficients 1 < ERS < 6 . In the lack of a better name, it is called the normal group. In this set, random fluctuations in GE space seem to play the main role in the genesis of cancer, as originally claimed in Ref. 27 . Notice that this group is conformed by very different tissues, from the medulloblastoma to the colorectal adenocarcinoma.
Group II, with five points in the figure, include cases in which genetic or viral causes enhance the rate µ ′ . The ERS index exhibits very high values in this set.
The abnormal values of ERS for the 7 tissues (12 points) contained in Group III could have an immunological origin. Indeed, our body uses physical barriers, humors and immune cells in order to protect the tissues against infections caused by pathogens, which are the most common attacks. The combined effects of these factors guarantees immunity. In tissues where one factor is predominant, the others could be somehow depressed. On the other hand, the protection against tumors, which come form inside, that is originate in tissue cells, is mainly the responsibility of immune cells. In other words, in tissues where the role of immune cells is depressed at the expense of increasing barriers or other components, the relative cancer risk, and correspondingly the ERS factor, is increased.
Barriers are known to play a basic role in the protection of germinal cells 39 and the brain 40 against infections. The cellular component of immunity in these tissues is, in some way, depressed with the purpose of avoiding inflammation events. The relatively high values of ERS could be explained in this way.
By contrast, the inclusion of the Medulloblastoma in the normal group is probably related to regional differences in blood-brain barrier permeability 41 .
With regard to bones, it is known that immunity relies strongly on defensins 42 , possibly with a depressed role of immune cells. On the other hand, the thyroid is known to have a close cross-talk with the immune system 43 .   www.nature.com/scientificreports/ It's dysregulation is the cause of immune disorders. One may speculate that a low cellular response is needed in order to prevent dysregulation of the thyroid. The extreme case in this group is gallbladder non-papillary adenocarcinoma, with an index ERS = 1300 , the understanding of which is a real challenge. However, one can speculate that the cellular response is also depressed in the gallbladder, because of the strong microbicide character of the bile 44 .
On the other hand, the relatively low value of ERS for the small intestine adenocarcinoma (eight times lower than the reference) can not have other explanation than overprotection by the cellular component of the immune system. Indeed, the small intestine is a possible entrance door for the microbiota living in the colon, and as such it requires special protection. The mean value of microbes/gm experiences a jump from 10 4 to 10 11 as we cross from the ileum to the cecum 45 . Barriers can not be reinforced because of the reduced dimensions. Thus, perhaps the Paneth cells 46 , Peyer's patches 47 , and other structures concentrated in the distal ileum are the responsible for this additional protection.
Finally, there is a group of 2 tissues exhibiting abnormally high values of the ERS index, presumably related to external factors. One example is lung adenocarcinoma, for which the concurrence of radioactive Radon and smoking produces a 90-fold increase of the slope.

Concluding remarks
In the present paper, the time evolution of microstates representing small portions of a tissue are described as Levy flights in gene expression space. The small amplitude Brownian component is characterized by a radius D √ t , much less than the distance between the normal and tumor regions, R =x 1 − R n − R t . Only sporadic