Intense terahertz pulses inhibit Ras signaling and other cancer-associated signaling pathways in human skin tissue models

Terahertz (THz) radiation has shown unique advantages in biomedical applications for novel diagnostic technologies due to the high sensitivity to molecular structure and chemical concentration. However, emerging evidence shows that intense pulses of THz radiation can induce significant non-thermal biological effects that must be characterized. In human skin exposed to intense THz pulses, relatively large responses characterized by differential gene expression profiles are observed. These data are analyzed by signaling pathway perturbation analysis to predict phenotypic endpoints and dysregulatory effects on cancer-related processes. The activities of several important pathways that drive the initiation, development, and progression of many human cancers are predicted to be suppressed, and this effect is intensity-dependent. Some affected pathways are targets for current and emerging anti-cancer therapies. In particular, the activity of the Ras signaling and Calcium signaling pathways is predicted to be significantly inhibited. These results indicate the possibility of an additional therapeutic mechanism of intense THz pulses, due to the potential for targeted suppression of pro-mitotic activity in diseased tissue.


Introduction
Terahertz (THz) radiation is a form of non-ionizing electromagnetic energy between 0.1 and 10 THz (3-300 cm −1 ). These frequencies overlap with natural oscillatory dynamics of rotational modes of gas molecules, low-frequency vibrations in molecular solids, and the stretching/twisting modes of hydrogen bond networks [1,2]. External THz coupling to these oscillations provides a sensitive modality to excite and probe natural dynamic behavior in many systems of interest.
Of particular interest is research and the application of THz radiation in biomedical fields [1,3]. The sensitivity to chemical structure and concentration has led to the development of diagnostic medical applications such as THz pulsed spectroscopy [4] and THz pulsed imaging [5,6], for which the contrast mechanism is predominantly differential water content, cell density, and protein concentration [7]. Overviews of diagnostic THz applications over recent years are reviewed in [1,[8][9][10][11]. The high sensitivity to molecular/chemical structure allows for quantitative dielectric characterization of the imaging subject with excellent contrast between diseased and healthy tissue states without the need for additional contrast agents. Recent studies characterizing the current diagnostic capability of THz imaging indicate comparable or improved diagnostic capabilities relative to other emerging or conventional modalities [6,[12][13][14][15][16].
Due to its non-ionizing nature, diagnostic use of THz imaging/spectroscopy is generally regarded as a non-destructive probe. However, the possible interaction mechanisms imply that EM energy at THz frequencies can couple to important molecular/cellular structures, which at high enough intensity may alter the structural dynamics to the point of affecting associated function. Experimentally, distinct THz absorption modes have been found in isolated DNA and protein samples [17,18], which motivates the hypothesis that coupling of high-intensity THz radiation to these important structures may dysregulate higher-level biological processes (cellular, tissue, organ, or organism levels) that rely on molecular-level function. Several studies have reported THz-induced phenotypic changes at all levels of biological organization: THz exposures have been shown to disassemble actin filaments [19], increase membrane permeability [20], non-thermally induce differential gene expression in human stem cells, keratinocytes, and skin models [21][22][23][24], affect cellular differentiation in mammalian stem cells [25,26], activate acute inflammatory responses at the cellular [22], tissue [23,24], and organism [26] level, and induce aneuploidy in human fibroblasts [27]. Highly intense pulses (∼MW cm −2 ) of THz radiation have been shown to activate several proteins of the DNA damage response, and induce phosphorylation of the H2AX histone in skin tissue models [28]. This is an indication of genotoxic stress often used as a correlative marker for double-strand breaks, a severe form of DNA damage that often results in cell death. As radiation-induced DNA damage is a known therapeutic mechanism for several genetic disorders including cancer, these data and observations led to the speculation that intense THz pulses may find application in a therapeutic capacity by achieving the same therapeutic result via wholly distinct interaction characteristics with biological systems, and warrants further investigation into the potential effects and clinical uses.
Among the relatively small number of studies conducted, there is a large degree of variance in reports of biological effects induced by THz exposures. This is largely due to highly heterogenous sets of exposure parameters, biological systems, phenotypic endpoints investigated, analysis methodologies, and reporting practices, as reviewed in [29][30][31].
In this paper, human skin tissue models are exposed to an extended train of picosecond-duration THz pulses of varying intensities, and the resulting changes to global gene expression are measured. Following the signal pathway perturbation analysis workflow outlined previously [24], signaling pathways that are likely to be dysregulated are identified, and the intensity-dependent magnitude of the predicted pathway dysregulation is presented. Of the identified pathways, special attention is paid to pro-mitotic and pro-inflammatory signaling pathways that are often dysregulated in human cancer. THz-affected genes that predominantly drive the predicted biological dysregulation are identified, and the THz-induced dysregulation to cancer-related processes in skin is reduced to a subset of only 88 genes. Novel cluster-based analyses of the pathway edge matrix that account for the highly complex cross-talk between distinct signaling pathways provide an ordered list of key THz-affected genes that predominantly drive dysregulation to cancer-related signaling processes. These provide a convenient roadmap for designing targeted investigative assays for precise mechanistic studies that are more practical than global transcriptomics (e.g. immunofluorescence protocols that allow detailed spatiotemporal characterization of target gene/protein dynamics in single cells), and represent the best candidates for understanding the biological mechanisms underlying THz-induced dysregulation of pro-mitotic/pro-inflammatory signaling in human skin. We hope these data provide the research community with a useful roadmap for motivating targeted experimental designs, or the discovery of novel therapeutic applications of intense THz pulses.

Signaling pathways in cancer
In cellular systems, signaling pathways are a series of biochemical interactions that serve to communicate changes to or within the cellular environment and regulate higher-level biological function. For example, extracellular ligands such as growth factors may bind to a receptor protein, initiating a serial cascade of protein binding/activation through the intracellular space that eventually activates transcription factors in the nucleus to carry out cell division. Abnormal activity of signaling pathways are associated with diseases such as cancer, and controlled modulation of signaling pathways are an effective method of treatment for many disorders. For example, the design of molecules that bind/inhibit proteins that participate in mitotic signaling pathways are popular chemotherapeutic agents. Signaling pathways can be represented visually by network graphs, in which each node represents a gene/protein, and the edges represent the empirically determined interaction characteristics between relevant proteins. This representation allows for prediction of the pathways and processes likely to be affected by an external agent (e.g. THz irradiation) from measurement of changes to the cellular system's gene expression profile. In this paper, several cancer-related pathways are identified as significantly dysregulated by intense THz pulses, based on the differential expression of gene transcripts that encode for the pathways' constitutive proteins. To explicate the above general description of signaling pathways, two pathways of particular relevance to cancer research and therapy are described, the Ras signaling and Calcium signaling pathways.

Ras and calcium signaling in cancer
The Ras signaling network is a ubiquitous and highly-conserved pathway in almost all eukaryotic cells, and regulates proliferation, survival, growth, migration, and differentiation [32][33][34]. The Ras protein, a membrane bound small GTPase of which there are three dominant forms (KRAS, NRAS, and HRAS), is a central node that regulates mitotic signal transduction by acting as a binary molecular switch to regulate Ras signaling activity. Due to this pathway's generalized mitotic function, it is one of the most commonly implicated pathways in the initiation, development, and progression of many human cancers [32]. Approximately 30% of all human cancers and 25% of melanomas are associated with a mutation in the central Ras protein [32,33]. These mutations allow Ras to acquire resistance to hydrolysis by GTPase-activating proteins (the switch's 'off ' signal), leading to perpetual pro-mitotic signaling that is a key feature in malignant transformation and metastatic invasion. For this reason, Ras and the related signaling components have been a key focus of cancer drug discovery research, with an emphasis on identifying and characterizing agents that inhibit the proliferative activity of Ras and the pathway's related signaling molecules [32][33][34].
Of additional interest is the THz effect on the Calcium signaling pathway, as calcium ions (Ca 2+ ) are one of the most potent second messenger particles in biological systems. In this pathway, Ca 2+ ions bind and activate calcium-binding proteins (CBPs) that regulate a diverse set of general cellular function such as protein signaling, metabolism, enzyme function, cytoskeletal dynamics, proliferation, and apoptosis [35,36]. In skin, epidermal differentiation is tightly regulated by a strong calcium gradient across the epidermis that triggers sequential expression of structural/binding proteins for stability and formation of the skin barrier. However, Ca 2+ ions may also bind/activate CBPs that drive tumor growth and metastasis through oncogene activation or dysregulated proliferation, apoptosis, adhesion, or morphology. Control of calcium signaling dynamics has therefore attracted significant interest as a potential target for novel cancer therapies [36,37].
As will be shown, intense THz pulses induce gene expressions in human skin tissue models that are predicted to suppress the signaling dynamics of these, and related, signaling pathways. In particular, the predicted suppression of pro-mitotic pathways that are often over-active in human cancer provides an additional avenue of potential therapeutic application of THz radiation. The therapeutic mechanisms proposed here involve THz-induced modulation of cell signaling in skin, and are distinct from the potential therapeutic mechanisms related to genotoxic stress and DNA damage proposed previously ( [23] and [28]).

Generation and detection of intense THz pulses
Detailed descriptions of the generation/detection methods for this work are outlined in [24]. Briefly, a 1 kHz train of single-cycle, picosecond-duration THz pulses was generated by optical rectification of infrared (λ 0 = 800 nm) laser pulses in LiNbO 3 [38], followed by black polyethylene to filter residual infrared/visible wavelengths. The pulse energy and spatial intensity distribution were measured with a pyroelectric detector (Spectrum Detector, Inc.) and camera (Electrophysics PV320), respectively. The temporal waveform was detected by free-space electro-optic sampling in GaP. The peak electric field was calculated by calibration to dielectric properties of the L = 200 µm GaP detection crystal (n 0 = 3.18, r 41 = 0.88 × 10 −12 m V −1 ) as where t GaP = 0.46 is the GaP transmission coefficient and ∆I/2I 0 is the measured peak modulation of the electro-optic signal as described in [24]. Equation (1) was used for the highest field calculation before the system was aligned for biological exposure, and lower field strengths were determined by scaling as the square root of the pulse energy, which has been verified to be accurate for this source. Five exposure conditions were used by attenuating THz transmission through crossed wire-grid polarizers, resulting in pulse energies of 0.03-2.4 µJ, peak electric fields of 27-240 kV cm −1 , and a measured spot size (1/e 2 ) of 1.6 × 2.7 mm 2 (additional details are provided in table 1). The THz pulse and spatial intensity distribution for the highest intensity pulse are shown in figure 1(a), and the corresponding frequency spectrum is shown in figure 1(b).

Preparation and exposure of skin tissue samples
Full-thickness 3D human skin tissue models were acquired from MatTek Corporation (EpidermFT, www.mattek.com/products/epidermft). These are a stratified co-culture of normal epidermal keratinocytes and dermal fibroblasts in a collagen matrix, and interact to form a mitotically and metabolically active 3D model of human skin. 3D tissue models have significant advantages over 2D cell monolayers as they allow the study to probe effects that arise due to more realistic intercellular interaction, cell-extracellular matrix interaction, and the impact of the micro-environment [39]. Upon delivery, tissue wells were transferred to  2.5 ml of pre-warmed media (DMEM, 10% FBS, 1X penicillin/streptomycin). The proprietary wells allow the tissue to absorb media from below through a porous membrane, leaving the epidermal layer exposed to air. This is particularly advantageous for THz exposure experiments, as this system allows for direct tissue exposure and bypasses constraints of attenuation losses when exposing through substrates or cell media. Tissues were incubated in this state for 16 h for re-equilibration. Tissue wells were removed from the incubator and placed vertically in a custom sample holder. To ensure precision of sample placement, the sample holder was first fit with a 1 mm alignment aperture and positioned such that THz energy is maximized, ensuring the center of the tissue well is repeatably coincident with the THz focus. The THz beam was focused to the epidermal layer for 10 min (figure 2(a)) to simulate realistic skin exposure conditions. The bulk/average temperature increase in the irradiated region on the surface of these samples as measured with a thermal imager (Reed Instruments, R2100) is less than 1 K, and these are benchmarked with similar measurements in a pure water reference environment. Following exposure, wells were transferred to fresh, pre-warmed media, and returned to the incubator for 30 min to allow for acute (early) response of THz-induced gene expression. To isolate the variation of THz pulse energy as the sole independent variable, similar exposure/fixing times as previous investigations were chosen [23,28]. The biochemical indicators of this response were fixed via snap-freeze in liquid nitrogen and transferred onto dry ice to −80 • C storage. This process was repeated for five different sets of THz pulse parameters (table 1) and one identical sham-exposed condition with the THz beam fully blocked. Each condition (5 exposed + 1 sham) was performed in quadruplicate for a total of 24 tissue samples.

Whole-genome expression profiling
The exposed region of the frozen tissue samples was excised with a 1.5 mm diameter sterile punch tool aligned to the center of the tissue wells to ensure only the THz-exposed cells were analyzed. Total RNA was isolated (Zymo Research, Irvine), and UV spectroscopy quantified RNA quality and integrity as per the manufacturers' instructions (NanoDrop, Wilmington, DE; Agilent 2100 Bioanalyzer, Santa Clara, CA). Amplified cDNA was labelled (Encore Biotin, IL) and hybridized to the HumanHT-12_v4 Whole Genome Expression BeadChip Arrays (Illumina). Differential gene expression and corresponding p-values were determined from the fluorescence measurements (iScan & BeadStudio, Illumina). P-values were corrected for multiple hypothesis testing via the false discovery rate method. Data was corrected for batch effects via the ComBat algorithm [40].

Signal pathway perturbation analysis
Pathway IDs, associated gene lists, and pathway topology information (i.e. interaction characteristics) were obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) [41]. Pathways were filtered for relevance to human cancer (KEGG categories include Environmental Information Processing, Cellular Processes, Organismal Systems, and Human Diseases). Signal pathway perturbation analysis was performed with the R/Bioconductor package ROntoTools, which allows for analysis of all datasets utilizing p-value-based weights on perturbation terms [42,43]. This algorithm combines conventional over-representation analysis with topology-based analyses that additionally account for the magnitude and direction of gene expression (upregulation vs. downregulation), as well as how genes interact within a given signaling network. This results in a prediction of the magnitude/direction of the perturbation status (activation vs. inhibition) of considered pathways, and provides information regarding potential phenotypic endpoints likely to be induced by intense THz pulses based on the gene expression measurements. Details on the theory and implementation of this algorithm can be found in [24]. Briefly, perturbations as measured by gene expression profiles are propagated through signaling networks that regulate known biological functions, taking into account the interactions (edges) between each gene (node/vertex). For each gene g i in the network, a perturbation factor PF (g i ) = ∆E (g i ) + A (g i ) is assigned, where ∆E (g i ) is the measured differential gene expression ratio (log-transformed), and the accumulated perturbation is estimated from upstream perturbation factors as where β ij describes the interaction direction between the i th and j th gene in the network (+1 for activation/anti-inhibition, −1 for inhibition/anti-activation, or 0 for no interaction), N ds is the number of genes downstream of g j , and the sum index of 'j < i' indicates that signaling perturbation may only come from upstream nodes. In this way, all network nodes are iteratively assigned a value of A (g i ). From this, a total accumulated perturbation score A tot = ∑ i A (g i ) is determined for the considered pathway, along with associated p-values adjusted for multiple hypothesis testing. Significantly perturbed pathways (p < 0.05) are considered dysregulated by the relative magnitude A tot . Positive (negative) values of A tot imply the associated pathway activity is increased (decreased) by intense THz pulses.

Identification of 'THz targets'
Not all genes differentially regulated by THz pulses (or any agent) affect a given pathway uniformly. For example, in a given pathway a node with many downstream effectors (i.e. a greater number of edges extending from the node) is more 'important' for signal proliferation, and induced pertubration A(g i ) of the node will contribute to total accumulated perturbation A tot moreso than genes with relatively few downstream effectors [42,43]. Moreover, a single gene may have distinct regulatory roles across multiple pathways. Hence, signaling pathways in nature are rarely separate closed systems, but rather there is a potentially large degree of cross-talk between sets of pathways that share common signaling proteins. It is therefore of interest to identify the genes that are both differentially expressed by THz and significantly drive pathway dysregulation, accounting for inter-pathway crosstalk between signaling systems containing common nodes. These represent genes that most significantly drive THz-induced changes of signaling processes in human skin, and are the best choice for determining the gene-level mechanisms that lead to the predicted THz-induced effects to cancer-related signaling. These provide a convenient roadmap to motivate future hypotheses and experiments in THz-cancer research. To achieve this, the 'pathway edge matrix' is constructed, for which the identified networks/pathways (nodes/genes) comprise the row-space (column-space). The matrix is populated by the normalized number of edges extending from the node in the corresponding network, and is a measure of the relative importance of that node (as nodes with a larger number of downstream effectors have greater potential to propagate node-level perturbation throughout the network). Standard k-means clustering (k = 2) is performed across rows/columns to group genes/pathways by matrix similarity. This algorithm iteratively bisects each dimension by maximizing the difference of group means. Since the metric in the matrix is a measure of gene significance in regulating cancer processes, the resulting gene list represents the best gene-level candidates for understanding mechanisms of the changes to cancer signaling dynamics induced by intense THz pulses.

Results
It is important to state that the measurements, analyses, and interpretations in this paper are relevant for the stated THz pulse parameters (figure 1, table 1) and biological system (human skin) under the stated exposure conditions (section 3.2). These may be generalizable to broader experimental parameters and cell-/tissue-types, but this must be verified rather than assumed.

Intensity-dependent global gene expression and signaling pathway dysregulation
To simulate a realistic skin exposure scenario, the THz beam was focused to the epidermal layer as shown schematically in figure 2(a). The nominal thickness of the epidermal layer is 65 µm, and is further comprised of a ∼30 µm layer of cornified keratinocytes ('corneocytes') and 35 µm active keratinocytes below. The effective penetration depth for the broadband THz pulse in water is calculated to be 75 µm. Using these, we estimate that the THz energy absorbed in the cornified epidermis, active epidermis, and dermis to be 33%, 25%, and 42%, respectively. Figure 2(b) shows the intensity-dependent global gene expression induced by exposure to intense THz pulses as described in the methods section. As THz pulse energy increases, the gene expression profiles display a general increase in expression magnitude and statistical significance. A relatively large genomic response was measured: At the highest pulse energy, 1681 genes met the conventionally defined significance criteria, with 1088 genes downregulated and 593 upregulated. From these gene expression profiles, pathway-level dysregulation was predicted as described in section 3.4. Of the 184 pathways available in the KEGG database, 54 were found to be significantly dysregulated by at least the highest THz pulse energy (p < 0.05). Of these, 8 are involved in regulating initiation, development, and progression of human cancer. These pathways are listed at left in figure 2(c), and the magnitude, direction, and significance of the pathway-level dysregulation are displayed by plots at center and right. These results indicate that intense THz pulses are predicted to dysregulate many important cancer-related signaling pathways. This predicted dysregulation generally depends on pulse energy, although the exact nature of the relationship remains to be elucidated. Figure 3 summarizes the results of identifying the most significant genes in THz-induced cancer signaling dysregulation, as described in section 3.5. As shown in figure 3(a), of the 1681 differentially expressed genes that meet significance thresholds, all dysregulation related to cancer signaling processes is due to a subset of only 88. These can be broadly categorized as in figure 3(b) as 'Pro-mitotic' , 'Cancer-specific' , and 'Pro-inflammatory' , for which the latter is seen to be nearly entirely genomically distinct. Figure 3(c) is the pathway edge matrix constructed as described in section 3.5, and shows the 88 genes within the eight cancer-related pathways predicted to be dysregulated. Since the matrix is populated by a measure of nodes' perturbative capacities, clustering these values leads to grouping of genes based on relative importance across pathways.

Transcript-level THz targets
An important feature of figure 3(c) is the structure that emerges in the clustered matrix that highlights the degree of inter-pathway crosstalk. The pathways are not independent signaling systems, but rather are regulated by many common genes in a complex web of inter-pathway crosstalk. This crosstalk is evident in the discrete blocks that emerge in the structure of the clustered matrix: Blocks along rows indicate groups of genes that significantly regulate a given pathway, and blocks along columns indicate individual genes that regulate function across multiple pathways. These groups are qualitatively indicated by the dendrograms at right (for pathways) and above (for genes). Due to these complex interactions, standard methods that identify dysregulatory sources within an individual pathway are insufficient, as the effect to other unconsidered pathways are not easily predictable but may dominate the higher-level phenotype. Therefore, the complete analysis accounting for complex inter-pathway crosstalk as described in section 3.5 is a statistical approach that allows one to identify primary sources of pathway-level dysregulation, taking all inter-pathway signaling contexts into account.
For row-clustering (pathways), there are three broad pathway categories that emerge as expected from the Venn diagram in figure 3(b): pro-mitotic, pro-inflammatory, and cancer-specific. The pathways belonging to each are given by the row-space of the edge matrix shown at left. These categories, the corresponding pathways, and predicted total accumulated perturbation at the highest THz intensity (A tot from figure 2(c)) are listed in table 2.
For column-clustering (genes), several discrete structures emerge in the clustered edge matrix as indicated by the pink box in figure 3(c). The ordering of the clustered column-space indicates the genes that are maximally responsible for the predicted THz-induced dysregulation to cancer signaling. Table 3 shows the first five gene clusters, as well as associated processes and THz-induced expression measurements at the highest THz intensity. In tables 2 and 3, bold (italic) text represents downregulation/inhibition (upregulation/activation) of the gene/pathway for the highest THz intensity.

Discussion
In earlier studies using similar skin models, THz exposure parameters, and biological assays as those discussed in this study, Titova et al observed differential expression of 442 genes induced by intense THz pulses [23]. Subsequent analysis showed a downregulation of genes in the epidermal differentiation complex and predicted activation of an acute inflammatory response. The present study verifies both observations in the relevant tissue system (human skin) and further characterizes specific pathways, associated biological processes, and the THz-affected genes that most significantly drive the predicted dysregulation to important cancer signaling processes.
As shown in figure 3(b), the predicted THz-induced dysregulation of 88 THz-affected genes in 8 cancerrelated pathways can be broadly categorized as activation of pro-inflammatory processes (positive A tot of the Cytokine-cytokine receptor interaction pathway) and suppression of pro-mitotic processes (negative A tot of the Glioma, Ras signaling, Rap1 signaling, Regulation of actin cytoskeleton, Calcium signaling, and cGMP-PKG signaling pathways). The suppression of pro-mitotic pathways, particularly within the Ras signaling and Calcium signaling networks, implies a potential therapeutic mechanism of intense THz pulses.

Activation of pro-inflammatory signaling
Cytokines are extracellular ligands that are released following inflammatory stimuli and participate in inter-cellular communication to regulate the growth, differentiation, and activation of immune cells [44]. Specific respones are triggered upon cytokines binding to cytokine receptors in the membranes of target cells. In addition to regulating the general inflammatory response in normal tissue, cytokines have an important role in the development of the tumor microenvironment and have been observed to both suppress and drive the metastatic phenotype, dependent on the cytokines involved or the type/stage of the cancer [44]. Emerging evidence shows that cytokines are involved in regulating tumor formation, and controlled modulation of this regulatory function may be exploited for cancer therapy [44,45]. Early phase clinical trials investigating multiple components of the immune system as therapeutic targets have shown success in effective tumor response in neural tissue as reviewed in [46].

The cytokine-cytokine receptor interaction pathway
The Cytokine-cytokine receptor interaction pathway is predicted to be the most strongly dysregulated by intense THz pulses, and is consistent with previous observations of a THz-induced inflammatory response [23,26,28,29,[47][48][49]. 18 of 76 genes considered in this pathway are differentially expressed at the highest THz pulse energy, and these are largely localized to differential expression of the chemokine and inflammatory cytokine families (CCL, CXCL), interleukins (IL6), and interferons (IL24), as shown in Cluster 5 of table 3. The predicted activation at the pathway level is due to accumulated perturbation at downstream nodes, which predominantly include the corresponding receptors of these ligand families (CCR, CXCR, ILR). Therefore, the dominant sources of pathway activation were not due to genes directly affected by THz exposures, but rather due to downstream targets of directly affected genes. This indicates that from a biological perspective, the extended exposure of skin tissue to intense THz pulses is recognized and responded to as an acute immune response. Kim et al have compared the tissue response of mouse skin to varying inflammatory stimuli and show that this inflammatory response to THz exposures is most genomically similar to a wound response, and dissimilar to responses stimulated by burns or exposure to other types of radiation (UV, neutrons) [26]. In this dataset, the activation of the inflammatory response is nearly entirely genomically distinct from the suppression of pro-mitotic signaling as shown in figures 3(b) and (c).

Suppression of pro-mitotic signaling
Although the Glioma pathway is predicted to be the most strongly suppressed in figure 2(c) and table 2 (A tot = −36.8, p = 0.004 at highest THz pulse energy), it is categorized as 'cancer-specific' , and the predicted dysregulation is due to the distinct effect on sub-networks, in particular the Ras signaling and Calcium signaling pathways. As evident in table 3, these are the most significantly dysregulated processes, due to downregulation of a relatively small number of genes (most significantly KRAS, MAPK3, and CBPs). The most strongly affected general cancer pathway is the Ras signaling network:

The Ras signaling pathway
The Ras signaling pathway, in particular the canonical RAS-RAF-MEK-ERK axis, is one of the most well-characterized signaling cascades in molecular biology, due to its ubiquity in nearly all eukaryotic cell types. As this pathway regulates general division, differentiation, and proliferation processes, it is often implicated in many human cancers and is therefore a target of modern cancer therapy research. The Ras protein, a membrane bound small GTPase of which there are three dominant forms (KRAS, NRAS, and HRAS), is a central node that regulates mitotic signal transduction by acting as a binary molecular switch to regulate Ras signaling activity.
Of the 100 genes considered in this network, 29 were significantly expressed at the highest THz pulse energy. THz-affected genes in the Ras superfamily include members of the Ras (KRAS, RAP1B, RRAS2), Rap (RAP1B), and Rho (RAC2) families. Associated Ras-dependent families (RASA1, RASA2, RASGRP1) additionally regulate the binding and activity of Ras-related proteins. As shown in table 3, expression changes of these genes (KRAS and MAPK3) are the most dominant source of cancer-related dysregulation and predicted suppression of mitotic signaling. Many of the gene clusters associated with the Ras signaling pathway are also responsible for the predicted suppression of several other networks in the pro-mitotic group (Rap1 signaling, Calcium signaling, cGMP-PKG signaling, and Regulation of actin cytoskeleton), indicating the nested nature of these signaling processes that interact with significant crosstalk. Processes that this group of pathways regulate, such as division, differentiation, motility, and apoptosis, are all involved in various aspects of mitotic activity. The THz-downregulation of the KRAS oncogene, associated signaling genes, and subsequent prediction of inhibition of Ras signaling activity (A tot = −26.1, p = 0.01), indicates that intense THz pulses may find potential therapeutic application for some cancers, with the goal of targeted inhibition of pro-mitotic signaling in diseased tissue.

The calcium signaling pathway
Eleven of 54 genes in this pathway were significantly expressed at the highest THz pulse energy, and these are largely represented by downregulation of calmodulin and calmodulin-like genes that encode for calcium binding proteins (CBPs; CALM, CALML, CAMK2D families, see Cluster 2 of table 3). These genes dominantly affected in the pro-mitotic category are directly related: many Ras-related proteins encoded by THz-affected genes are in turn calcium-regulated, and their proper function relies on regulatory activity of the listed CBPs [50]. Since calcium influx may be expected to dysregulate division/differentiation pathways that are closely correlated to calcium signaling, inhibition of associated signaling may be a genomic response to field-induced membrane permeation as discussed in [20]. As differential gene expression profiles represent the tissue response to the applied stimulus, expressing genes to attenuate calcium signaling is a feasible homeostasis-seeking response to increased calcium influx. This hypothesis is additionally corroborated by the pathway analysis in figure 3(c), since the pro-mitotic pathways identified are also largely calcium regulated. For example, RAS activity is sensitive to changes in Ca 2+ concentration, since several Ras-activators (e.g. Rap1) are in turn calcium-activated via binding of RASGRP [51]. In addition, studies investigating gene expression in skin following electroporation have reported activation of a local inflammatory response, including increased expression of the gene families upregulated by intense THz pulses (CCL and CXCL families) as discussed in section 4.1 [52].

Consideration of thermal effects
An important consideration is the possibility of thermal effects in biological systems induced by THz radiation, as discussed in [29,31,47,[53][54][55]. While recent studies have shown that thermal effects induced by relatively low powers of CW THz exposure is genomically distinct from general bulk heating processes [56], care has been taken in these experiments to ensure that the observed biological response is non-thermal by limiting the pulse train's duty cycle (1 kHz train of picosecond-duration pulses). This minimizes the average power of the THz beam such that thermal effects are negligible, while maintaining high peak power for THz energy to propagate in high-attenuation aqueous environments [53].
Two effects are occurring: (a) each pulse will cause some amount of temperature increase in a volume of space, and (b) the thermal energy is dissipated from this volume at a rate dependent on the sample's thermal conductivity. These two contributions for the 1 kHz pulse train will reach a steady state of average heating. An estimate of the per-pulse temperature increase is ∼4 mK (using ∆Q = mc∆T, where ∆T is the temperature change due to a THz pulse with total energy ∆Q in a volume of water of mass m = ρV (ρ = 1 g cm −3 ). Volume V was calculated using measured THz spot size and penetration depth and c = 4.2 J g −1 K −1 is the specific heat capacity of water), and the average steady-state heating due to the 1 kHz pulse train was measured with a thermal imager in water to be less than 1 • C. This measurement is consistent with heating measurements induced in similar THz exposure systems when scaled by average power [21,54,56]. In terms of biological response, the heat shock factor (HSF, regulator of the heat shock response) was not differentially expressed, and other heat shock proteins were not upregulated. No biological processes related to thermal stress/response were significantly identified by pathway or gene ontology analyses. It is therefore unlikely that the observed gene expression induced by intense THz pulses is thermal in origin. At minimum, intense THz pulses were not recognized biologically as a normal thermal stimulus.
An additional possibility that should be explored in future work is dependence on the pulse train's repetition rate (1 kHz for this study). As argued above, both the per-pulse and train-averaged temperature increase are biologically negligible. However, a small thermal energy input, while generally biologically insignificant, may induce a frequency-dependent effect by coupling to and altering structural/chemical dynamics. In addition to intense THz wave interaction with fast biomolecular motions, additional thermal kHz-coupling to slower dynamics may be another key interaction mechanism to explore as well.

Conclusion
In this paper we show that intense THz pulses cause significant changes in gene expression patterns in human skin. The expression magnitude and statistical significance generally increase with THz pulse intensity. From these data, intensity-dependent dysregulation of cancer-related signaling pathways are predicted, and eight cancer-related pathways are significantly identified at at least the highest THz pulse energy. Pro-inflammatory processes are predicted to be activated due to upregulation of common inflammatory cytokines, interleukins, and interferons, and is consistent with previous studies' observations. Pro-mitotic pathways are predicted to be inhibited, largely due to downregulation of genes in the Ras superfamily or genes that encode for calcium binding proteins. Pathway topology information was used to identify 88 genes (of 1681 possible candidates) responsible for all predicted cancer-related dysregulation, with 42 genes dominantly responsible. These are presented as an ordered list in table 3, along with the associated signaling pathways, and represent the strongest candidates for future specific mechanistic studies in characterizing THz-induced dysregulation of the signaling processes of interest. Suppression of the Ras signaling and Calcium signaling pathways are of particular significance, as these are commonly implicated pathways in human cancers, and popular targets for existing and emerging cancer therapies.
These observations suggest a potential therapeutic mechanism of intense THz pulses. Conventional cancer therapies involve inhibition of proteins that are responsible for uncontrolled cell division/proliferation (e.g. by pharmacological agents specially designed for the target oncoprotein) or driving metastatic cells towards the apoptotic phenotype (e.g. by irreperably damaging DNA with bleomycin or ionizing radiation). However, many of these therapies induce significant side effects, such as the absorption of ionizing radiation or systemic activity of anti-mitotic drugs in otherwise healthy tissue. Due to high absorption and very low scatter, non-ionizing THz exposures have a very high degree of conformality to the target region. Additionally, since applications are limited to areas of direct radiation delivery (skin, oral, endoscopic, intra-operative, etc.), exposure to normal tissue is minimal. THz exposures may therefore be less susceptible to certain disadvantages of current therapies. More work detailing observations in many different sample types, in particular focusing on the differential effect in normal cells relative to cancer cells, must be done to isolate these potential effects that may be exploited clinically. These include a systematic variation of relevant THz exposure parameters such as dose, spectral content, or repetition rate as discussed above, and standard viability assays to investigate prolonged effects in daughter cells, or to relate the predicted signaling dysregulation to changes in the cellular phenotype. Microarray-based approaches as presented here are resource-intensive and only provide a broad snapshot in time of the global genomic response. This study provides an empirically-obtained ordered list of THz-affected genes that dominantly drive signaling dysregulation, and we hope this provides the research community with valuable guidance in choosing targets and forming precisely-designed assays (e.g. utilizing standard protocols such as immunofluorescence that characterizes spatiotemporal expression dynamics with high detail). These methods are additionally much less intensive on lab resources, and may be performed at high-throughput levels necessary for characterizing spatiotemporal signaling dynamics of target genes at the functional protein level. If THz-regulation of key cancer signaling genes show a significant, predictable response that is therapeutically effective, it may offer a viable option for some cancers in the future.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.