Levy model of cancer

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, and reveals an interesting connection between the risk and the way the tissue is protected against infections.

Each point in this figure comes from a biopsy, Small samples are taken off from different patients and processed in order to obtain expression values for 60483 genes. For each gene, we define a reference value, e ref , by geometric averaging over the normal samples. Then, new variables are defined:ê = log 2 (e/e ref ). The origin of coordinates in this figure is precisely the center of the cloud of normal samples,ê = 0. A covariance matrix is defined and diagonalized. The first eigenvectors are used to define new coordinate axes: PC1, PC2, etc. Details may be found in Ref. [2]. We shall only stress that the first axis, PC1, which accounts for 51 % of the total data variance, is the cancer axis, which allows the discrimination between normal and tumor states.
Let v 1 be the eigenvector of the covariance matrix along PC1, andê the expression vector corresponding to a given sample. Then, x 1 =ê · v 1 is the position along PC1 of the sample. Normal samples define a region around the origin with r.m.s. radius R n = 11.71. On the other hand, the cloud of tumor samples is centered atx 1 = 155.89, and its r.m.s. radius is R t = 28.53 [6].
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 [6,7]. 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. We have computed in Ref. [7] 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 [8]. The conclusion is that the tumor minimum is the deepest in Fig.  1, the one with highest fitness. Our drawing for the fitness distribution is a sketch, but we are convinced that at least a rough quantitative description can be reached.
The intermediate region, R n < x 1 <x 1 − R t , holds a low-fitness barrier [6,7], 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 [6].
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 as 1.5 × 10 7 [9]. 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 [10]. 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. 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 [11]. Second, there is also a rate of accumulation of epigenetic (mainly methylation and phosphorylation) events modifying the normal expression levels [12]. Both processes could be boosted by inherited mutations [13,14] or external carcinogens [15].
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 gene i. Eq.
(1) describes a Markov chain of events [16]. 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 [2]. v 1i > 0 and v 1i < 0 correspond, respectively, to genes that should be over-and under-expressed (silenced) in order to the tumor to progress. We have distinguished the genes CST1 and AQP8. The former is a known marker of colon cancer [17], whereas the latter plays a significant role in colon homeostasis [18] 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 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 they cancel. In this way, Eq. (1) for the small displacements in GE space describes a 1D Brownian or Poisson process [19].
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 Err 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 Err(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 proto-tumor could be the subject of an attack by the immune system in the very early stages [20]. By definition, the constant is less than zero because the overall constant in Eq. (5) is less than one.
In Table 1 we compile a set of parameters for a group of tumors. The geometry of the normal and tumor regions, i.e. the parametersx 1 , R n and R t come from Ref. [6]. The D value is estimated as the maximum of |v 1i | [2]. 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. [21,22]. The reported values of risk represent averages over 380 cancer registries from different cities and countries in the world [22].
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 1. 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. [21] 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 × 10 −6 ), 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 small amplitude random variations in GE values. In the next section, we shall consider large jumps.
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 [23,24]. However, also as a consequence of aging there is an accumulation of epigenetic events and DNA damages leading to a reduction of fitness and a displacement towards the low-fitness zone. Thus, aging acts in the same direction as the low amplitude fluctuations of GE values.
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 [25].
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 [26]. The global states of these networks define attractors. 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 mutations and GE variations leading to cancer, triggered by let unknown causes, which is the basic hypothesis in the atavistic theory of cancer [27].
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 [28] or Levy [29] kind, with a power-like tail. Indeed, the Pareto character of GE distribution functions was demonstrated in Ref. [30] (see also [6]). The Levy character of the length distribution functions in mutations was shown in [31].
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 [31], Levy flights are known to take place in many other biological processes, for example foraging [32].
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 = 115.65 with the scale D = 0.0526. 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 due to mutations 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 COAD GE distribution functions the exponents take values between 1.6 and 2.0 [6]. 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 1: The constant should be negative according to our hypothesis of µ small. The results of the test are shown in Fig. 4. The observed behavior is consistent with a linear dependence with slope near one. The Pearson correlation coefficient is 0.85, with a p-value equal to 0.04. The main reason for the remaining dispersion of points could be the assumption that the rate of large jumps, µ, is roughly the same for all of the tumors. A 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 Eq. (11) in order to re-examine the data presented in paper [21]. The qualitative idea is to set a reference value for the coefficient in front of the r.h.s. of Eq. (11), a ref = 2 × 10 −14 , and include an extra risk score (ERS), in such a way that: We should try to understand the observed values of ERS in terms of the tissue characteristics. The results are shown in Fig. 5. In order to facilitate the analysis, the studied tumors are separated in groups.
Group I includes 11 tumors (9 tissues), located in a band delimited by red dashed lines in Fig. 5, and coefficients 1 < ERS < 5. 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. [21]. 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 8 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.
Barriers are known to play a basic role in the protection of germinal cells [33] and the brain [34] 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 [35].
With regard to bones, it is known that immunity relies strongly on defensins [36], 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 [37]. 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 [38].
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 [39]. Barriers can not be reinforced because of the reduced dimensions. Thus, perhaps the Paneth cells [40], Peyer's patches [41], and other structures concentrated in the distal ileum are the responsible for this additional protection.
Finally, there is a group of 3 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 mi-crostates 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 large jumps, of Levy nature, allow the microstate to reach the cancer basin of attraction.
Although it is understood that aging induces a motion in the direction of the low-fitness region, it was not explicitly included in our model. Work along this direction is necessary.
The resulting formula for the risk of cancer in a tissue was quantitatively tested against observed data, and applied to the qualitative analysis of the risk of cancer in a set of tissues. The most important conclusion, in our opinion, is a possible connection between the risk and the way the tissue is protected against infections. The blood-brain barrier in the cerebrum, for example, preventing the entrance of pathogens, is also the reason for the relatively low rate of elimination of prototumors, and thus large risk per stem cell in this organ. The low risk per stem cell in the small intestine, on the other hand, is understood as a reinforcement of the cellular component of immunity.     Table 1 is used to this end. A nearly flat dependence on the abcisa is obtained, thus small amplitude fluctuations in gene expression space may not account for the risk of cancer in these tissues..