A multiscale model of the role of microenvironmental factors in cell segregation and heterogeneity in breast cancer development

Abstract

We analyzed a quantitative multiscale model that describes the epigenetic dynamics during the growth and evolution of an avascular tumor. A gene regulatory network (GRN) formed by a set of ten genes that are believed to play an important role in breast cancer development was kinetically coupled to the microenvironmental agents: glucose, estrogens, and oxygen. The dynamics of spontaneous mutations was described by a Yule-Furry master equation whose solution represents the probability that a given cell in the tissue undergoes a certain number of mutations at a given time. We assumed that the mutation rate is modified by a spatial gradient of nutrients. The tumor mass was simulated by means of cellular automata supplemented with a set of reaction diffusion equations that described the transport of microenvironmental agents. By analyzing the epigenetic state space described by the GRN dynamics, we found three attractors that were identified with cellular epigenetic states: normal, precancer and cancer. For two-dimensional (2D) and three-dimensional (3D) tumors we calculated the spatial distribution of the following quantities: (i) number of mutations, (ii) mutation of each gene and, (iii) phenotypes. Using estrogen as the principal microenvironmental agent that regulates cell proliferation process, we obtained tumor shapes for different values of estrogen consumption and supply rates. It was found that he majority of mutations occurred in cells that were located close to the 2D tumor perimeter or close to the 3D tumor surface. Also, it was found that the occurrence of different phenotypes in the tumor are controlled by estrogen concentration levels since they can change the individual cell threshold and gene expression levels. All results were consistently observed for 2D and 3D tumors.

Fig 1. Gene regulatory network for breast cancer.

Microenvironmental agents, estrogen (pink) and oxygen (blue) affect gene states transcription. The arrows indicate activation whereas the short bars represent inhibition interactions between genes.

Fig 2. The structure of the multiscale model for breast cancer.

The diagram illustrates the processes that take place at each spatial scale. The colored arrows illustrate the coupling between the scales. At each spatial scale the corresponding equations that were solved are indicated.

Fig 3. Gene expression levels in terms of the normalized concentrations of oxygen (O) and estrogen (E) for three initial values of the activation/inhibition (A-I) threshold parameters.

(A) ${\eta }_{0_i}=0.25$ low, (B) ${\eta }_{0_i}=0.5$ medium, and (C) ${\eta }_{0_i}=0.75$ high.

Fig 4. Microarray showing gene expression and cell phenotypes.

Phenotypes were identified with attractors in the state space of the GRN dynamics. Normal cell phenotypes are represented by the states a and c while precancer phenotypes are represented by the states b, d, e, f, g, h, and i. Cancer phenotype was represented by the state j, which corresponds to the over-expression of genes HER2, ATK1, P21, and CDK2. The dotted line of the state h means that there is an overlap with state i. The initial threshold parameters for each of the ten genes from top to bottom are: η_oi = (0.3, 0.45, 0.15, 0.45, 0.15, 0.45, 0.3, 0.15, 0.3, 0.5). Phenotypes were identified through the analysis of the GRN boolean dynamics as well as from an analysis of the dynamics of a GRN continuous model. (For details see Section B in S1 Appendix).

Fig 5

(A) Spatial distribution of cell phenotypes and (B) segregation index for a 2D tumor. These results were obtained for parameters values: P(B) = 1, κ_ϵ = 0.1, ξ = 0.5, χ = 4, α₁ = α₂ = α₃ = 8 × 10⁻³, λ₁ = 100, λ₂ = 50, and λ₃ = 200.

Fig 6. Spatial distribution of the number of mutations for (A) 2D and (B) 3D tumors.

These results were obtained for the same parameters values as in Fig 5.

Fig 7. The microarrays show the fraction of cancer cells in a 2D tumor for (A) Estrogens expression levels versus genetic inheritance values.

(B) Estrogens expression levels versus epigenetic inheritance values. (C) Estrogens expression levels versus estrogen receptors probability. (D) Genetic inheritance versus epigenetic inheritance. The microarrays were obtained with the following parameter values: (A) P(B) = 1 and χ = 0; (B) P(B) = 1 and ξ = 0; (C) χ = 4 and ξ = 0.5; (D) P(B) = 1 and κ_ϵ = 0.1. The consumption parameter values were α₁ = α₂ = α₃ = 8 × 10⁻³, λ₁ = 100, λ₂ = 50, and λ₃ = 200.

Fig 8

The microarrays corresponding to a 2D tumor and represent: (A) Mean segregation index. (B) Fraction of normal cells in a tumor, (C) Average number of mutations. (D) Shannon index that measures tumor heterogeneity. The axes show low to high gradient levels of the genetic and epigenetic contributions. They were obtained for the following parameter values: P(B) = 1, κ_ϵ = 0.1, α₁ = α₂ = α₃ = 8 × 10⁻³, λ₁ = 100, λ₂ = 50 and λ₃ = 200.

Fig 9

Averages carried out over various system realizations for different quantities as a function of time evolution for several values of the genetic and epigenetic inheritances parameter values that represent moderate epigenetic (Moderate EP) and high epigenetic levels (High EP) as well as low genetic (Low G) and moderate genetic levels (Moderate G). The remaining parameter values are the same as in Fig 8. (A) Percentage of precancer cells, (B) Percentage of cancer cells, (C) segregation index mean value, and (D) activation-inhibition threshold mean value for gene expression. Thin lines correspond to averages over 2,4,6…20 system realizations. Thick lines represent the averages over 10 system realizations.

Fig 10. Microarrays corresponding to a 2D tumor.

(A) Fraction of cancer cells in the tumor. (B) Tumor shape and spatial distribution of segregation index. (C) Fractal index, D) Shannon index. They were obtained for the following parameter values: P(B) = 1, κ_ϵ = 0.1, ξ = 0.5 and χ = 4.

Fig 11

Microarrays obtained for 3D tumors that show: (A) Fraction of cancer cells in the tumor, (B) The spatial distribution of mutations for different values of α₃ and λ₃. (C) Average threshold values and (D) tumor heterogeneity. The axes show low to moderate gradient levels of the genetic and epigenetic contributions. The results shown in (A), (C) and (D) were obtained for the parameter values: P(B) = 1, κ_ϵ = 0.1, α₁ = α₂ = α₃ = 8 × 10⁻³, λ₁ = 100, λ₂ = 50, and λ₃ = 200.

Fig 12. Spatial distribution of phenotypes.

(A) 2D tumor and (B) 3D tumor for different estrogen consumption rate values α₃ and supply rate values λ₃. (C) 2D tumors representation for different values of oxygen consumption rate values and (D) 2D tumors for different glucose consumption rate values. In all figures we set the model parameters as in Fig 11, except in panel (C), where we set α₁ = α₃ = 16 × 10⁻³ and λ₁ = λ₃ = 100. In panel (D), we set α₂ = α₃ = 16 × 10⁻³ and λ₂ = λ₃ = 100.