Ergodicity Breaking and Self-Destruction of Cancer Cells by Induced Genome Chaos

1. Introduction

The concept of ergodicity plays a central role in many problems of dynamical systems. Ergodicity is defined as the equivalence of time averages and phase space (or ensemble) averages in a physical system. Besides its theoretical significance, ergodicity is crucial for practical aspects of measuring the system’s parameters, since sufficiently long observations of temporal behavior are often impossible and, therefore, it is desirable to substitute them with finite time measurements [1,2,3]. In biophysical and biomedical applications, ergodicity is particularly important for processes that may involve chaotic behavior. These processes include the transport of blood and other liquids, intra- and extracellular fluid flow in cytoplasm and nucleoplasm, and multiphase flow involving macromolecular biopolymer solutions [3,4,5,6].

Typically, the chaotic behavior implies that a dynamic system is ergodic, since the information about the system’s history is constantly erased due to growing fluctuations. However, the interrelation between chaos, turbulence, and ergodicity is often quite complex. Experimental simulations suggested that turbulent chaotic behavior governed by the Navier-Stokes equation is ergodic [7]. Despite that, many biological systems are non-ergodic. The causes of ergodicity breaking are diverse. Many of them simply lead to the anomalous diffusion deviating from the classical linear Einstein–von Smoluchowski dependency of the mean-square displacement on the lag time. Such causes include macromolecular crowding, aging-induced flow through obstacles, and the so-called “hydrodynamic memory” [8]. In all these cases, the chaotic behavior is responsible for ergodicity breaking and transport deceleration in comparison with the classical diffusion law.

Ergodicity breaking is also related to the creation of information and entropy production in a dynamical system. In an ergodic system information (or memory about past states) is constantly erased, and therefore entropy is created.

While ergodicity breaking in dynamical systems is associated with the chaotic trajectories in their configurational space, the concept of genome chaos has been suggested to explain rapid and uncontrollable evolution of malignant cancer cells [9,10,11,12]. Simple models with three competing cell populations (host, immune, and tumor cells) may result in complex chaotic behaviors [13,14]. Chaotic behavior can lead to a significantly higher maximum tumor size when compared to non-chaotic behaviors, since increasing the parameter associated with the killing of tumor cells by immune cells is demonstrated to increase the maximum tumor size [15,16]. Nikolov et al. [17] found that the growth and progression of tumors has features similar to “strange attractors” in the dynamical systems combining local instabilities with global stability. Rocco et al. [18] suggested that the emergence of distinct growth phenotypes in clonal populations is related to weak ergodicity breaking.

Heng et al. [12] proposed that karyotype-mediated macroevolution, rather than gene mutation, is the common driving force for most cancers. Their model implies a two-phased cancer evolution mechanism with genome alteration-mediated macroevolution followed by gene mutation/epigenetic alteration-mediated microevolution. High stress induces genome chaos and macroevolution allows for new system survival, while microevolution aids system modification, providing an advantage for stepwise proliferation and competition [12].

One particular type of genome chaos is associated with chromothripsis [19], a massive genomic rearrangement, which is of particular interest for the formation of cancer. Chromothripsis is a catastrophic event when the chromosome is shattered into slices and then reassembled in a random order with mutation rates significantly increased (Figure 1a). Tossing of chromosomal fragments is associated with ergodicity.

In the present paper, we will suggest an analytical model for entropy creation and ergodicity during chromothripsis, we will also present a computational model for the mutation of the BCL2 gene responsible for apoptosis in cancer cells. In addition, we will review experimental data on time and population averages of cloned organotropic supermalignant RA-2 rhabdomyosarcoma of rats stem cells, and investigate the possibility and experimental evidence that increasing the degree of chaos beyond the non-ergodic threshold may lead to the self-destruction of the tumor cells.

2. Materials and Methods

In this section, we present models for calculating entropy and ergodicity defect in the presence of chromathripsis as well as the rate of mutation of the BCL2 gene responsible for cancer cells’ apoptosis.

2.1. Model for genome chaos

In physics, chaos is defined as a property of a complex system whose behavior is unpredictable due to its great sensitivity to small changes in initial conditions. Usually, such systems are dynamically unstable, and small fluctuations grow exponentially with time. While the concept of genome chaos is widely used, there is no common definition of this term. A process of complex, rapid genome re-organization, caused by chromosomal instability and resulting in the formation of chaotic, unpredictable genomes, is usually implied by genome chaos [10,11,12].

We apply the concepts of ergodicity and entropy to genome chaos. First, we assume that the mutation rate μ(θ), defined as the probability of mutation per nucleotide per division depends on a quantitative mutagenicity parameter, θ, characterizing the pressure from the environment increasing the probability of mutations.

During chromothripsis, up to thousands of chromosomal rearrangements can occur in a single event (Figure 1a). Thus, if the original composition of an affected chromosome is $\left(G_1,G_2,\dots G_N\right)$, where G_n is the n-th fragment of the chromosome, after shattering and consequent stitching, the sequence of genes can rearrange in a random manner as $\left(G_{m_1},G_{m_2},\dots G_{m_N}\right)$. Some genomic information can be lost during the event. According to estimates, DNA breaks repaired by mitotic gene conversion are accompanied by surprisingly high mutation rates which are more than 1000-fold higher than spontaneous mutations [19]. The rate of spontaneous mutations in somatic human cells is about ${\mu }_0=1.4\times {10}^{-10}$nucleotides per cell per division [20]. In cancer cells without chromothripsis such rate may be on the order of ${\mu }_{can}={10}^{-7}$ [21].

The exact mechanisms of shattering and stitching remain obscure; however, there is growing evidence that the decrease of protein P53, which plays a central role in maintaining genome stability, correlates with chromothripsis. The corresponding gene, TP53, is the most frequently mutated gene in human cancer [22]. The appearance of chromosomal bridges and micronuclei in cells serves a visual manifestation of chromothripsis (Figure 1b). As far as the frequency of chromothripsis in cancer cells, in a recent survey of 4,934 cancers, Zack et al. [23] suggested that chromothripsis occurred in 5% of all samples, with frequencies ranging from 0% in head and neck squamous carcinoma to a maximum of 16% in glioblastoma.

We now will distinguish three stages of the behavior of the system. At the first stage, no chromothripsis occurs. Following Rocco et al. [18], the cell state is viewed as a point in a multidimensional configuration space, characterized by gene expression profiles, $x_n$, where the epigenetic landscape is a hypersurface formed by inverse probability for a cell to be at a certain state. The epigenetic landscape is analogous to the energy landscape of a Hamiltonian system, $H\left(x_n\right)$. A phenotype is defined as the basin of attraction of each stable state (local minimum) on the landscape. The average barrier height between the stable states, ∆h, characterizes the probability of switching between the states (in other words, of a mutation) and can be estimated as $∆h=1/{\mu }_0=7\times {10}^9$. The average time the system spends at a certain phenotype’s basin of attraction, t₀, is larger than the observation time, for that reason, the system is essentially non-ergodic [18]. Essentially, the cells are in homeostasis.

In the second stage, massive rearrangements occur due to chromothripsis. Rearrangement of the fragments results in the mixing, which can be estimated quantitatively using Shannon entropy [24]

$$S_r\left(N\right)=-N\mathrm{l}\mathrm{n}\left(N\right)$$

(1)

For N chromosomal segments participating in the rearrangement, the number of rearrangement variants is P(N)=N! A simple statistical physics analogy can be a phase transition caused by increasing temperature above the melting point. The stability in a thermodynamic system at constant temperature and pressure is characterized by the Gibbs free energy defined as the difference between the enthalpic and entropic terms. The positive sign of the Gibbs energy change prohibits the spontaneous reaction

$$∆G=∆h-\theta ∆S>0$$

(2)

where the enthalpic term, ∆h, characterizes the resistance to mutation. Hence, ${\theta }_0=\frac{∆h}{∆S}$. Note that the purpose of Eq. 2 is to establish an analogy with the thermodynamics of phase transitions rather than to determine the exact value of the mutation rate at which the chromothripsis occurs. This is because the quantitative characteristics would depend on the definition of the variable θ.

During chromothripsis, essential mixing of the chromosomal segments occurs. The degree of ergodicity can be estimated by calculating the ergodicity defect. For ideal mixing, the process is expected to be ergodic and the state of the cancer cells is not homeostatic.

At this stage, the mutation rate exceeds the critical threshold, $\mu >{\mu }_{cr}$, is above the critical threshold ${\mu }_{cr}$, so that mutations are too intense for cells to survive. As far as the biochemical mechanism, the gene that is responsible for apoptosis in cancer cells is BCL2 [25]. It is therefore hypothesized that when the critical number of mutations in cancer cells is achieved, apoptosis is induced.

Several measures of deviation from the ergodic behavior have been introduced in the literature to account for the non-ergodic behavior. Földes-Papp and Baumann [5] (2011) suggested decoupling the effects of the molecular crowding and the temporal heterogeneity by presenting the power exponent, which controls the dynamics of the interaction network, as a product of these two factors. Scott et al. [26] suggested the ergodicity defect D, defined at different scales (on a map T) with respect to a basis of functions f given by an integral of the square of space and time averages

$$D\left(f,T\right)\propto \int {\left(f^{\mathrm{*}}\left(x,T\right)-\overline{f}\right)}^{2}dx$$

(3)

where $f^{*}$ and $\overline{f}$are the time and space averages.

2.2. Gene mutation and translation

The BCL2 cDNA sequence (ENST00000398117.1) was retrieved from the COSMIC database. A Python script was utilized to simulate genetic mutations with the following parameters: a mutation rate of 10^-4 and 60 cell divisions (Hayflick limit). The gene sequence was initialized as a list of characters, each representing a specific DNA sequence for the simulation. The apply_mutation function was used to randomly alter bases in the gene sequence according to the mutation rate. For each base, there was a probability of mutation, where the base could change to any nucleotide (A, T, C, or G). A genetic code dictionary was employed to translate both the mutated gene sequence and a reference gene sequence into protein sequences. A translation function, translate_dna_to_protein, utilized this dictionary to convert both the reference and mutated gene sequences into protein sequences.

2.3. Protein stability prediction

The wild-type (wt) structure, a component of the Bcl2-BINDI complex, was sourced from the Protein Data Bank as a crystal structure (PDB ID: 5JSN) with a resolution of 2.1 Å [27]. Predictions for mutant forms of the Bcl-2 protein were generated using the SWISS-MODEL algorithm [28] according to the protocol published elsewhere [29]. To assess protein stability, the standard Rosetta protocol was implemented, involving the calculation of the energy score (ES). This stability prediction method was executed by aligning root-mean-square-deviation (RMSD) values with Rosetta energy parameters, as outlined by Ramelot et al. [30]. The differences in Rosetta energy scores for mutated forms (ΔES) were computed using the following equation:

$$\Delta ES={ES}_{mut}-{ES}_{wt}$$

(4)

where ES_mut and ES_wt are the Rosetta energy score values for the wild type and mutated forms, respectively.

2.4. Protein Function Prediction

The crystal structure of Bcl-2 in association with a Bax BH3 peptide was retrieved from the Protein Data Bank (PDB ID: 2XA0) with a resolution of 2.7 Å [31] to determine the Bcl2-Bax binding site. The molecular interaction between the protein and the peptide was evaluated utilizing the ZDOCK molecular docking server for the prediction of protein-protein complexes and symmetric multimers [32]. The binding energy (E_bind) was calculated by using the PPI-affinity tool designed to predict and optimize the binding affinity of protein-peptide and protein-protein complexes [33].

3. Results

In this section, we present modeling results for mutation of the BCL2 gene as a main cancer pro-apoptotic inhibitor during chromothripsis and experimental observations on increasing genome chaos by selection [34,35,36].

3.1. Modeling results for the rate of mutation of the BCL2 gene

Before the stability and function of Bcl-2 were analyzed, the wild-type BCL2 gene (comprising 720 nucleotides) underwent a random mutation, resulting in a gene sequence with a Hamming distance of 1. This mutated sequence was subsequently translated into a protein sequence, which exhibited a sequence identity of 99.58% and a similarity score of 1277 when compared to the wild-type protein. A single nucleotide polymorphism (SNP) at amino acid position 5 was detected, with a G→W substitution.

Figure 2. Energy Score (ES) in Rosetta Energy Units (REU) vs. RMSD scatter (a,b) and contour plots (c,d) for the protein molecules including wild-type and mutants of Bcl-2. The energy thresholds (ES = -280 REU and RMSD = 0.5 Å) are depicted as dashed lines.

Figure 2. Energy Score (ES) in Rosetta Energy Units (REU) vs. RMSD scatter (a,b) and contour plots (c,d) for the protein molecules including wild-type and mutants of Bcl-2. The energy thresholds (ES = -280 REU and RMSD = 0.5 Å) are depicted as dashed lines.

Following this, a Rosetta protocol based on the Monte-Carlo method identified a significant decrease in stability (ΔES = 4.11 REU) for the mutated protein. The energy score (ES_mut) for the G5W mutation was -278 REU, compared to -282.11 REU (ES_wt) for the wild-type protein (refer to Figure 1b,c). The energy score for the mutated protein exceeded the energy threshold (TS = -280 REU), indicating the effect of the SNP. All minimum-energy conformations were within the threshold range (0.5 Å) determined in our previous experiments for the molecular folding of the Trp-Cage protein [37,38]. Figure 3 presents the 3D structural alignments of the wild-type and mutant template-based homology models with the RMSD value of 0.46 Å, the wild-type homology model alongside its minimal structure with the RMSD value of 0.31 Å, and the mutant homology model and its minimal structure with the RMSD value of 0.27 Å.

Figure 3. The 3D structural alignments of the wild-type and mutant template-based homology models (a), the wild-type homology model alongside its minimal structure determined by the Monte-Carlo method (b), and the mutant homology model with its minimal structure determined by the Monte-Carlo method (c).

Figure 3. The 3D structural alignments of the wild-type and mutant template-based homology models (a), the wild-type homology model alongside its minimal structure determined by the Monte-Carlo method (b), and the mutant homology model with its minimal structure determined by the Monte-Carlo method (c).

To evaluate the mutation’s impact on protein folding, the kernel density estimate was employed to compare energy score versus RMSD plots, which often form a triangular shape indicative of a protein folding pattern. The base of the triangle (lower RMSD values) represents near-native structures that have been refined and thus have lower energy scores. Consequently, the wild-type protein exhibited a well-formed triangle, while the mutated form displayed a distorted folding pattern, suggesting the mutation’s influence on protein folding (Figure 2c,d). Finally, the interaction between Bcl-2 and a Bax BH3 peptide was assessed to determine the impact of the SNP on the inhibition of apoptosis in the cancer cell. The results indicated that the Bcl-2-Bax interaction (Figure 4A,B) was impaired by the point mutation due to an increase in the binding energy for the mutated form (E_bind = -11.1 kal/mol) compared to the wild-type (E_bind = -11.3 kcal/mol).

3.2. Case study: self-destruction of rhabdomyosarcoma (RA) of rats

Already in the 1990s, Kravtsov and co-workers studied transplantable rhabdomyosarcoma (RA) from rats induced by 20-methylcholanthrene in lung tissue. They conducted several stages of selection for increasing genome chaos [34,35,36]. The frequency of cells with micronuclei (FCM) was used as a marker of genome chaos. Experimental metastases-clones RA-2 had a spherical shape with dimensions of 3 - 5 mm. Micronuclei were usually observed in single-nucleus cells in interphase, as evidenced by their external morphological features. The size of the observed micro-nuclei varied from 1 to 3 μm (Figure 5).

Artificial selection of cell clones with a high degree of spontaneous karyotype instability was used. Clones 3–5 mm in diameter prepared from the lungs of grafted rats were cut in half, one half of each clone was used to prepare a smear, and the other half of the clone was placed in a vial with medium 199. As a rule, the duration of the determination of FCM for a sample of 50–60 clones did not exceed 6 hours. The selection for increasing FCM was effective and led to a gradual increase in FCM to 6.9 - 8.1% (Table 1). After 3-4 steps of selection for an increase in FCM, more than 25% of the clones had an FCM greater than 9.0%. The average value of the FCM was 6.9%, and the range of variability was from 0.5 to 17.0%.

Metastases-clones stopped their proliferation and were subjected to apoptosis at the average FCM close to 4.7% and the maximum FCM was close to 15%, roughly speaking, on the order of 10%.

4. Discussion

Chromothripsis increases the rate of mutation in cancer cells by three orders of magnitude leading to the rates of ${\mu }_c~{10}^{-4}$mutations per division per nucleotide. If in the population of n cancer cells, the fraction of cells with chromothripsis constitutes f_c, the average mutation rate is given by

$$\mu =f_c{\mu }_c+\left(1-f_c\right){\mu }_{can}$$

(4)

Assuming the values of $f_c=0.1$ corresponding to 10% FCM, ${\mu }_c~{10}^{-4}$ and ${\mu }_{can}~{10}^{-7}$ yields $\mu =1.009\times {10}^{-5}$ mutations per division per nucleotide.

The pre-selection distribution of chromathripsis cells (Minimum FCM = 0 %, Average FCM = 0.6%, maximum FCM = 2.0 %) is assumed to be an equilibrium distribution. The population of cells maintains this distribution throughout the generations of division. Therefore, the time-average (generation-average) value is equal to the population average and the population is ergodic. In other words, the clone-average value does not change with time and therefore coincides with the time-average. However, during selection of cells for increasing chromothripsis the ergodic equilibrium is shifted away from the ergodic state. Non-ergodic state eventually leads to the destruction of all cancer cells by apoptosis.

The effective mutation rate was calculated with Eq. 3 and the ergodicity defect was calculated using Eq. 4. The values are presented in Table 1 and Figure 6.

The data suggest that during the selection process, the number of cells affected by chromathripsis grows quickly reaching maximum possible values at the fifth stage, or even earlier in some experiments [34,35,36]. Increasing the frequency of cells with micronuclei (i.e., likely subject to chromathripsis) results in a sharp increase of the ergodicity defect. The process continues until the concentration of FCM reaches 10%-15% of total cancel cells.

Modeling the BCL2 gene mutation for 60 generations (Hayflick limit) at $\mu ={10}^{-5}$ mutations per division per nucleotide yields the overall probability of mutation at probability on the order of p=1.

We, therefore, conclude that the threshold value of the concentration of cells with chromothripsis corresponding to the deterioration of Bcl2 or a similar protein responsible for genomic stability is on the order of $f_c=0.1$ (or 10%). This may explain the observation that the RA-2 cells tend to reach the maximum concentration of cells with chromothripsis before being subjected to apoptosis. The observed phenomenon provides potential new ways of cancer treatment by inducing tumor cells’ apoptosis.

The transition to genome chaos leads to ergodic mixing of genome segments, whereas the further increase of the mutation rate results in apoptosis of the cancer cells (Figure 7). To some extent, a physical analogy can be suggested with the transition from solid to liquid state of matter (with a chaotic motion of molecules) and then to the evaporation with increasing temperature. Increasing mutation rate plays the role of increasing temperature in this case.

5. Conclusions

Genome chaos is associated with cancer propagation. Chromothripsis is one of the most prominent manifestations of genome chaos. During chromothripsis, mixing of the chromosomal segments occurs and, therefore, entropy is created while information is erased. The regime with erasing information, i.e., without memory, is ergodic. However, with the increasing rate of mutations, malignant cells deviate from the ergodic regime and, eventually, apoptosis happens. Computational modeling of the mutations in the BCL2 gene responsible for the apoptosis in tumor cells shows that the probability of mutations approaches the unity at the concentration of the cells with chromothripsis on the order of 10%. This is consistent with experimental data on cloned rhabdomyosarcoma from rat stem cells, artificially selected for increasing rate of genome chaos. Clones with the frequency of chromothripsis above 10% result in cell apoptosis possibly due to mutations in the BCL2 gene. Potentially this can be used for suppressing cancer cells by stimulating increased mutation rates and shifting them into a non-ergodic proliferation regime.