Introduction

As applications of terahertz (THz) radiation rapidly emerge1,2,3,4,5 researchers are beginning to explore the effects of THz radiation on biological matter and it is becoming increasingly evident that THz radiation can influence cellular processes in various ways6,7,8,9,10,11. Recently, it was experimentally demonstrated that THz radiation can influence gene expression and accelerate stem cell differentiation toward adipocytes12,13. A whole variety of effects, including changes in gene expression, mitotic disturbances14 and DNA production15 were detected in Chinese hamster ovaries15, human-hamster hybrid cells16, eukaryotic unicellular organism17, rat18, mouse12,13 and human19 cell cultures as well as, very recently, in E. coli cultures20 after treatment with millimeter wave non-ionizing radiation. Apart from revealing the potentially detrimental effect of THz radiation, a detailed investigation of the effects of THz radiation as an external stimulus, capable of influencing cellular gene expression, could be of great practical importance. The ability to selectively and remotely activate specific genes would be a revolutionary change with tremendous technological and economical implications in biotech industries that rely on protein or lipid production. Moreover, with control of gene expression comes control of cellular function, structure and cell differentiation (reprogramming), which will revolutionize regenerative medicine and offer great medical advances21.

In this study, building on our previous theoretical and experimental work12,13,22,23, we investigate the specificity of THz irradiation on mouse mesenchymal stem cells (mMSCs). In particular, we used i) an ultra-short pulsed broadband (centered at ~10 THz) tabletop source24,25 to irradiate mMSC cultures for 2 and 12 hours and ii) a single-frequency (SF), 2.52 THz, continuous wave laser source to irradiate mMSC cultures for 2 hours. In the experiments, we controlled and measured the temperature and found that the temperature-increases were not significant. Correspondingly, the expression levels of genes encoding heat shock proteins were practically unaffected.

Our previous findings12,13 revealed that prolonged (9-hour-duration) broadband THz irradiation of mMSCs applied 120 hours immediately after having treated the cells with differentiation medium inducing adipose phenotype26, can accelerate cell differentiation by activating the transcriptional cascade regulator PPARG27. In the present experiments, we chose to irradiate cells for 2 and 12 hours with different THz sources at an earlier differentiation time point (48 hours after treatment with the differentiation medium) in order to assess the robustness of the above-mentioned effect as well as evaluate the specificity of gene expression due to irradiation parameters. Similar to previous findings12,13, the irradiation conducted at this new (earlier) time point results in overexpression of genes closely related to the adipose phenotype. However, in contrast to our prior work12,13, no morphological changes were observed. Interestingly, while our previous experiments did not reveal differentially expressed genes after only a 2-hour-duration of THz irradiation (with either broadband or SF), here our microarray survey and RT-PCR experiments show that there are indeed some differentially expressed genes, albeit different than the ones observed after 9 or 12 hours of THz irradiation. In particular, we found that genes affected by prolonged irradiation are characteristic for already differentiated cells, i.e. for adipocytes, while genes differentially expressed after short (2 hours) THz irradiation are primarily characteristic of pluripotent stem cells. Further, our microarray analysis implies that irradiation with SF laser or broadband THz source affects different genes. These findings suggest that the THz effect on gene expression depends on the duration and type of THz radiation, as well as on the stage of stem cells differentiation prior to irradiation. Finally, our DNA dynamics simulations demonstrate that there is strong DNA-breathing activity in the core promoter regions of differentially upregulated genes identified by our microarray survey. This association is consistent with our previously proposed mechanism of THz-DNA interaction22,23.

Results

We performed three different irradiation scenarios. Two of the scenarios involved using a broad-spectrum (centered at ~10 THz) pulsed tabletop THz source25 applied to irradiate mMSCs for 2 hours and 12 hours, respectively (Fig. 1a). In the third scenario, a single frequency, 2.52 THz, continuous wave laser source was used to irradiate mMSCs for 2 hours. In each experiment, a control mMSC culture was placed adjacent to the irradiated sample (screened from the THz radiation) and the temperature of the cultures was kept at 26–27°C for both the control and irradiated sample (Fig. 1b). Immediately after completing the irradiation, total RNA was extracted from the irradiated sample and the control and microarrays were used to identify differential changes in gene expression between the irradiated sample and its respective control. Each scenario was applied in three independent biological experiments and the reported results were averaged. In each experiment, the mMSC cultures were synchronized to be at the same differentiation time point immediately before the irradiation (Materials and Methods).

Figure 1
figure 1

The experimental setup for THz irradiation.

(a) Schematic depiction of the two-color generated broadband ultrafast broadband THz source including its spectrum with and without a Petri dish. The single-frequency continuous wave THz laser source (Edinburgh Instruments FIRL 100) is not shown but its frequency (2.52 THz) is depicted in purple on the spectrum in the inset. (b) Schematic depiction of the experimental setup displaying the irradiated and control mouse mesenchymal stem cultures. Both the irradiated sample and its control were held at the same temperature.

THz radiation drives mouse MSCs toward differentiation by affecting gene expression

First, we confirmed our earlier observations12,13 for each of the irradiation scenarios. Prolonged (12 hours) broad-spectrum THz irradiation of mMSCs resulted in overexpression of PPARG, adiponectin, GLUT4 and FABP4 (Fig. 2a), as previously observed12,13. In contrast, this prolonged irradiation, applied only 48 hours after the treatment with differentiation medium, did not cause the clear lipid droplet-like inclusions in the cellular cytoplasm that we have previously observed12,13 in mMSCs exposed to the same radiation conditions but irradiated at a later differentiation time point (120 hours after chemical treatment). This phenotypic difference suggests that the THz influence depends on the level of stem cell differentiation. Furthermore, irradiation of mMSCs for only 2 hours (regardless of the THz source) did not change the expression of PPARG, adiponectin, GLUT4 and FABP4 (data not shown), consistent with the previous reports12,13. Similarly, there was no elevation of gene expression of heat shock proteins in the irradiated samples (Fig. 2b). These results demonstrate an overall reproducibility in gene expression changes due to non-thermal mechanisms of THz radiation, while elucidating a new parameter that influences the THz effect – namely, the level of mMSC differentiation.

Figure 2
figure 2

Expression of adipose related, heat shock and stress related genes.

(a) RT-PCR measurements for selected transcripts from cells irradiated for 12 hours with a broadband THz source. The results confirm our previous observations12,13 that PPARG, adiponectin, GLUT4 and FABP4 were overexpressed in samples after prolonged irradiation, while PPARA remains unaffected. (b) RT-PCR measurements for selected transcripts of genes encoding heat shock and stress-related proteins from cells irradiated for 12 hours with a broadband THz source showed no overexpression in the irradiated samples.

Microarray analysis of the total RNA extracted from the mMSCs in each of the irradiation scenarios revealed statistical evidence (B-value > 0) for 20 differentially expressed genes (Table 1). In the cells irradiated for 12 hours with the broad-spectrum THz source 4 genes were overexpressed and 16 were underexpressed (Table 1). Importantly, Gem and Slco4a1 were underexpressed while Nfe2l2 was overexpressed. In addition to PPARG, adiponectin, GLUT4 and FABP4, these three differentially expressed genes provide further evidence for an accelerated differentiation of mMSCs when irradiated for 12 hours with a broad-spectrum THz source. Indeed, it was previously shown that Gem is expressed in the G1 phase of mouse pluripotent stem cell, while already differentiated cells lack Gem due to proteosomal destruction mediated by the anaphase-promoting complex/cyclosome (APC/C)28. Similarly, the membrane protein Slco4a1 is overexpressed in pluripotent stem cells29 but suppressed in differentiated cells and this protein has been used for cell sorting in flow cytometry experiments, as a marker for human embryonic stem cells30. In contrast to Gem and Slco4a1, which are suppressed in differentiated cells, Nfe2l2 is connected with adipogenesis31 and is overexpressed in adipose tissue32. Taken together, the observed underexpression of Gem and Slco4a1, combined with the upregulation of Nfe2l2, is suggestive of a higher level of differentiation in mMSCs irradiated for 12-hours with the broadband THz source as compared to (non-irradiated) control cell cultures.

Table 1 Differentially expressed genes. The list of genes identified by the microarray analysis as possessing statistically significant differential expression (B-value > 0) in each irradiation scenario

Furthermore, hierarchical clustering based on the 20 differentially expressed genes identified by the microarray analysis results in a clear separation between the irradiated and control samples (Fig. 3a). While there is minor heterogeneity between the biological triplicates, the overall changes of gene expression due to THz irradiation are consistent in all triplicates, thus allowing the formation of the two distinct groups. Interestingly, even in this small pool of differentially expressed genes, a gene ontology (GO) analysis reveals statistically significant (p-value < 0.001) overrepresentation of underexpressed genes in 2 GO terms and statistically significant (p-value < 0.001) overrepresentation of overexpressed genes in 5 GO terms (Fig. 3b). The 5 GO terms of upregulated genes provide an indication of a positive effect of THz radiation on the cell cycle machinery, while the 2 GO terms of downregulated genes hint at a lack of DNA damage stimulus and cellular stress.

Figure 3
figure 3

Analysis of differentially expressed genes in cells irradiated for 12 hours with broadband THz source.

(a) Heatmap for 20 differentially expressed genes (with statistical evidence B-value > 0) displaying the formation of two groups – one containing irradiated samples and the other containing their respective controls. Color depicts differential expression (control vs. irradiated samples) where red color corresponds to underexpression and bright yellow to overexpression. (b) Gene ontology analysis revealing overrepresentation of over- or underexpressed genes in GO terms for differential expressed genes identified in samples irradiate for 12 hours with a broadband THz source.

The microarray survey reveals only a few statistically significant differentially expressed genes after 2 hours of THz irradiation. Only three genes were affected by the 2-hour broad-spectrum irradiation: Slco4a1 was upregulated, while Car3 and Ear11 were downregulated (Table 1). Slco4a1 was also the only statistically significant differentially expressed gene (upregulated) in the cells irradiated for 2 hours by the single frequency THz source. Interestingly, Slco4a1 is overexpressed after 2 hours of THz irradiation (regardless of the THz source) while being underexpressed after 12 hours of broadband irradiation. One possible explanation of Slco4a1's expression patterns is that the high propensity for DNA breathing of Slco4a1's core promoter (see below) allows a catalytic effect of THz irradiation on its expression, when Slco4a1 is transcriptionally active (i.e., highly expressed), which occurs while the cells are still pluripotent; this effect diminishes as the mMSCs become differentiated (under prolonged THz irradiation), i.e., when other pathways become active and Slco4a1 becomes suppressed in the differentiated cells29.

Terahertz radiation has a parameter specific effect on gene expression

Our microarray survey has revealed a pool of genes with strong statistical evidence for differential expression in our three irradiation scenarios. Previous studies have demonstrated that some of these affected genes, namely Car3, Ctfg, Fn1, Nfe2l2, Gem and Slco4a, are involved in cell differentiation. For example, Car3 participates in PPARG regulation33 while Gem, Nfe2l2 and Slco4a are known markers for pluripotency. Furthermore, Ctfg has been implicated in fibroblast differentiation of human mesenchymal stem cells34, whereas Fn1 (a known pluripotency factor) is overexpressed in adipose-derived stem cells35. An RT-PCR examination of these genes is in agreement with our microarray analysis, while it also illuminates an inter-sample heterogeneity of gene expression (Fig. 4). Car3 and Nfe2l2 are underexpressed in cells irradiated for 2 hours (regardless of THz source) but overexpressed in cells irradiated for 12 hours with a broadband THz source. In contrast, Gem and Slco4a are upregulated in cells irradiated for 2 hours (regardless of THz source) while being downregulated in cells irradiated for 12 hours with the broadband THz source. Interestingly, Ctfg is overexpressed in cells irradiated with a broadband THz source (regardless of the duration of THz irradiation) while being underexpressed in cells irradiated with a single frequency THz source. Thus, Ctfg displays an opposite regulation when irradiated with different THz sources. Taken together, these inter-sample gene expression heterogeneities indicate that the effect of THz radiation (at least for the examined genes) depends on THz parameters such as the duration and type of radiation. However, such gene expression heterogeneity is not universal, as there are genes for which THz radiation, regardless of duration and laser source, has no effect: examples are genes encoding heat shock proteins (Fig. 2b).

Figure 4
figure 4

Comparison of the differential expression between genes involved in stem cell differentiation.

RT-PCR measurements (irradiated samples normalized to their respective controls) for selected set of genes in the three irradiation scenarios. The y-axis depicts the average gene expression fold change between controls and irradiated samples; as such a value of 1 corresponds to the same average expression in both controls and irradiated samples.

The RT-PCR experiments also indicate a limited sensitivity of our gene chip survey. While the microarray analysis correctly identified the differential expression for genes where strong statistical evidence exists (B-value > 0), it was unable to find all differentially expressed genes. For example, the microarray analysis accurately identified Gem as underexpressed after 12 hours of broadband exposure but it was unable to find that Gem was overexpressed in the other two irradiation scenarios, which was demonstrated by RT-PCR. This limited sensitivity in detecting differentially expressed genes could potentially be improved by using more replicates of irradiation experiments within restricted irradiated conditions (e.g., with narrow frequency intervals, smaller amplitudes of the THz field, etc.) followed by either microarray or RNA-seq analysis. Examining genes for which there is less statistically significant evidence (p-value < 0.05, not adjusted for multiple hypothesis testing) for their differential expression may provide a hint toward the homo- or heterogenic effect of the THz radiation on a larger scale. We found that these differentially expressed genes with p-value < 0.05 strongly differ in each of the irradiation scenarios, as shown in a Venn diagram (Fig. 5a). Interestingly, hierarchical clustering based on the differential expression of each of these three different sets of genes displays again formation of two groups – one of the irradiated and another of non-irradiated mMSCs, for each irradiation scenario (Fig. 5). Thus, in all three cases, one of the groups is formed by the irradiated triplicate mMSC cultures, while the other group consists only of their respective controls. This separation of the patterns of gene expression into two groups supports the idea of gene and radiation parameter specificity of the THz effect. Hence, although mMSC cultures have similar response to one type of THz radiation, the changes in duration or THz source dramatically affect differentially expressed genes. While only more irradiation experiments (followed by gene expression analysis) will unequivocally reveal the complete cellular heterogeneity of gene expression patterns due to THz irradiation, our results provide an indication that the effect of THz on gene expression has an intrinsic global heterogeneity, not limited to few genes.

Figure 5
figure 5

Global gene expression patterns of heterogeneity due to THz irradiation.

(a) Venn diagram of upregulated (green) and downregulated (red) genes with p-value < 0.05 (not adjusted for multiple hypothesis testing) in the three irradiation scenarios. Note that only one gene is shared between all three irradiation scenarios. (b) Heatmap for the differentially expressed genes (p-value < 0.05) identified in mMSCs irradiated for 2 hours with SF THz source. (c) Heatmap for all differentially expressed genes (p-value < 0.05) identified in mMSCs irradiated for 2 hours with the broadband THz source. (d) Heatmap for all differentially expressed genes (p-value < 0.05) identified in mMSCs irradiated for 12 hours with the broadband THz source. The heatmaps were generated using genes with p-value < 0.05 (not adjusted for multiple hypothesis testing) and each heatmap displays the formation of two groups – one containing irradiated samples and another containing their respective controls. Color depicts differential expression (control vs. irradiated samples) where red color corresponds to underexpression and bright yellow to overexpression.

Propensity for DNA breathing and gene expression changes caused by THz irradiation

While our experimental data indicates the existence of a THz radiation parameter specific effect on gene expression, the precise mechanisms governing this influence remain unclear. Our previous mesoscopic modeling of DNA breathing dynamics in a THz field22 suggest that, even at very low amplitude of the radiation field23, specific irradiation may affect gene expression by inducing local unzipping of DNA – i.e., by creating new open states in the double helix that in turn can affect transcription initiation36,37,38 or transcription factors binding39,40. Consistent with this modeling, we have previously found that genes with high propensity for breathing are overexpressed in response to THz radiation (e.g., the cascade transcription regulator PPARG27) while genes with low propensity for breathing remain unaffected by the THz radiation12 (e.g., PPARA). It is likely that the THz effect on gene expression is only catalytic and depends on many factors that are interplaying in the living cell, e.g., availability of transcription factors (suppressors or enhancers), position of the dynamically active promoters in the chromatin structure, as well as on the nature of the pathways that are active and/or activated by the THz irradiation, etc.

Here, we simulate DNA breathing dynamics of the core promoters of two upregulated pluripotency markers, namely, Gem and Slco4a1. Our evaluation of the propensity for breathing dynamics reveals the existence of large transient bubbles at the transcription start sites in the core promoters of each of these genes (Fig. 6). The dynamical activity of these promoters is consistent with our hypothesis that THz radiation can influence (in a nonlinear manner) the low frequency biomolecular vibrations of the intramolecular hydrogen bonds and drive new large-amplitude intramolecular motions that can change local stability of the double helix and in turn affect the transcription of their respective genes.

Figure 6
figure 6

Computer simulations of the propensities for breathing dynamics of the core promoters of two pluripotency markers.

EPBD-based Langevin molecular dynamics showing the intrinsic breathing dynamics in the core promoter regions of Gem and Slco4a1. The x-axis shows the base pair position in the core promoter with transcription start site (TSS) depicted as +1. The y-axis shows the bubble length (in base pairs) while the color corresponds to the average bubble probability for bubbles with amplitude bigger than 3.5 [Å].

Discussion

Our results show that THz irradiation of mMSCs can cause specific catalytic changes in cellular function that are closely related to the gene expression and differentiation state. The strictly controlled experimental environment indicates minimal temperature changes and the absence of any discernable response to heat shock and cellular stress, thus implying a non-thermal cellular response to the applied low power THz stimulus. Our RT-PCR experiments and gene expression survey suggest that some genes in the irradiated mMSCs cultures are activated, while others are suppressed. Importantly, the different irradiation conditions result in dissimilar regulation of Car3, Ctfg, Gem, Nfe2l2 and Slco4a, indicating that gene expression is dependent on the parameters and duration of THz irradiation. Interestingly, even 2 hours of irradiation are sufficient to cause differential changes in gene expression for certain active genes (e.g., Slco4a1).

In our previous work12,13, we irradiated mMSC cultures 120 hours after cell treatment with medium that induces differentiation toward adipose phenotype and observed activated expression of marker genes for adipocytes as well as lipid droplet-like inclusions in the cellular cytoplasm. Here, we report changes in gene expression of the adipocyte markers, but lack of morphological changes in mMSCs that were irradiated (at the same THz conditions) only 48 hours after the treatment with the same differentiation medium (i.e., presumably at an earlier differentiation state). Hence, irradiating stem cells further in their differentiation program (viz., ones irradiated 120 hours after treatment) results in noticeable morphological changes that are not displayed in cells irradiated at an earlier differentiation stage. This observation suggests that the effect of THz irradiation depends on the differentiation state of the stem cells.

One significant concern surrounding these results is that no comprehensive model can explain the effect of THz radiation. Our suggestion22,23 is that THz radiation may affect gene expression by interacting (in a nonlinear resonant way) with the hydrogen bonds of biomacromolecules, whose vibrational frequencies are in the THz range41,42, e.g., by perturbing the conformational dynamics of double-stranded DNA or by affecting protein-DNA binding. This suggestion is based on prior work37,43,44 that establishes a strong relationship between breathing dynamics of double-stranded DNA and cellular function36,38, as well as between DNA breathing and protein-DNA binding in some cases39,40,45. Our suggestion is also supported by the fact that millimeter wave radiation are capable of melting DNA oligomers by disrupting only intramolecular hydrogen bonds46. Furthermore, our experimental results are consistent with this proposed explanation, since some overexpressed genes have high propensity for DNA breathing. Nevertheless, the observed results may be due to other effects such as changes in the cell's membrane potential47, or protein-protein interactions48, or in the protein conformation itself49. Further molecular level experimental investigations of THz radiation's ability to induce specific openings of double-stranded DNA and changes in protein-DNA binding are needed to fully determine the validity of our hypothesis.

Finally, the performed microarray survey and RT-PCR experiments reveal a small group of differentially expressed genes in all three experimental irradiation scenarios. Even this small group of affected genes shows opposing effects of THz radiation on gene expression and points toward the need for comprehensively cataloging genes that can be affected under specific THz irradiation parameters. The generation of such a catalogue is not a simple task and it will involve a large-scale experimental effort with potentially tens of replicates irradiated under a large variety of conditions, taking into account (at least) parameters such as: amplitude of the THz field, narrow frequency intervals and exposure duration. Nevertheless, the generation of such a catalogue will be instrumental in opening the door to using THz radiation for non-contact control of cellular gene expression. We are positive that additional efforts, such as including protein and functional readouts as well as in vivo and ex vivo investigations of the effect of THz radiation are needed to validate, generalize and obtain stronger insight into the effect of the millimeter wave radiation on gene expression. Further insight into how the THz radiation can alter: the DNA breathing; or cell membrane potential; protein-DNA or protein-protein binding; or in general conformations of biomacromolecules (e.g., the local bending), is needed to shed light on the mechanisms of the THz effect on biological matter.

Methods

Mouse stem cell culture

Mouse mesenchymal stem cells (MSCs) (ScienCell Research Laboratories, CA 92011, Catalog Number M7500) were cultured on tissue culture treated plastic T-75 flasks (Corning). Once the cells reached 80–90% confluence, the cells were sub-cultured in supplemented medium (95% α-MEM, 5% fetal bovine serum (FBS) and antibiotics) for THz irradiation treatment. Roziglitazone (1 μM), insulin (5 μg/ml), 3-Isobutyl-1-methylxanthine (100 μM) and dexamethasone (1 μM), were added to the medium 48 hours prior to irradiation.

Terahertz sources and irradiation

Through the frequency mixing of ultrafast fundamental and second-harmonic laser fields at 800 nm and 400 nm, respectively, a directional electrical current in pressurized atomic gas (argon) plasma can be generated via two-color optical field ionization when the optical beams come to a focus inside the gas cell. In the case of ultrafast lasers (<100 fs) this technique is capable of producing electromagnetic radiation at THz frequencies50. The collinearly propagating second harmonic field is produced by placing a β-barium borate (BBO) nonlinear crystal inside the pressured gas cell. We use a recently developed coherent-control scheme optimized for such type of THz generation in gases, yielding a new source of high-energy (~1 μJ, pulse width 35 fs duration, i.e., high peak power per pulse ~30 MW), broadband THz radiation (centered at ~10 THz, λ = 30 μm, Fig. 1a) at a high repetition rate of 1 kHz25,50. In addition, we used a far-IR optically pumped molecular gas continuous wave THz laser source (Edinburgh Instruments FIRL 100). All experiments with this laser were conducted at frequency 2.52 THz (λ = 118.8 μm).

Analysis of affymetrix gene chip assay and RT-qPCR

Total RNA was individually extracted from the 9 irradiated cell cultures and their 9 respective controls and separately hybridized on 18 Affymetrix gene chips (GeneChip Mouse Genome 430 2.0 Array). The Affymetrix gene chips hybridization procedures were conducted in KUGR: the Genomic Center at the University of New Mexico, as recommended by the supplier (Affymetrix). The raw microarray data (.CEL files) discussed in this publication have been deposited in NCBI's Gene Expression Omnibus51 and are accessible through GEO Series accession number GSE41106.

The raw microarray .CEL files were analyzed using R and Bioconductor52. More specifically, the robust multi-array (RMA) method included in the affy package53 was used for preprocessing, background correction and normalization of the raw data. Probes without Entrez identifiers or gene symbols as well as probes with small standard deviation across samples were filtered out from the analysis. Differentially expressed genes were identified using the Bioconductor limma54 package and ranked by calculating the B statistics for each gene. Gene annotation was performed using the Bioconductor mouse4302.db package, while the hypergeometric test of the Bioconductor's GOstats55 package was used for identifying overrepresentation of differentially expressed genes in gene ontology terms.

Cellular RNA was extracted using the RNeasy Mini kit (Qiagen) following the manufacturer's instructions. For cellular qRT–PCR, 500 ng total RNA was reverse transcribed using the reverse transcribed (Clontech, Advantage RT-for-PCR kit) with random-oligo primers; qRT-PCR reactions were assembled in triplicates with 1 μM gene specific primers and 5 μl Applied Biosystems SYBR Green PCR Master Mix in 10 μl reactions. qPCR was executed in a 384-well block Applied Biosystem 7900-HT Real-time PCR machine. Gene specific primers were selected (PRIMER BLAST, NIH). PPARA and 18S were used as internal control. The comparative Ct method (Applied Biosystem)56 was used to analyze the data resulting from the RT-PCR experiments. All RT-PCR data are represented as mean ± standard error.

Computer simulations

Our extended Peyrard-Bishop-Dauxois (EPBD) model is an extension of the classical Peyrard-Bishop-Dauxois nonlinear model57 to include inhomogeneous stacking potentials36. Our EPBD-based Langevin molecular dynamics and Markov Chain Monte Carlo computer simulations and the analysis of the dynamic trajectories and DNA equilibrium states, from which we derived the average bubble probabilities, were performed as previously reported38,58.