Disassembly of microtubules by intense terahertz pulses

: The biological effects of terahertz (THz) radiation have been observed across multiple levels of biological organization, however the sub-cellular mechanisms underlying the phenotypic changes remain to be elucidated. Filamentous protein complexes such as microtubules are essential cytoskeletal structures that regulate diverse biological functions, and these may be an important target for THz interactions underlying THz-induced effects observed at the cellular or tissue level. Here, we show disassembly of microtubules within minutes of exposure to extended trains of intense, picosecond-duration THz pulses. Further, the rate of disassembly depends on THz intensity and spectral content. As inhibition of microtubule dynamics is a mechanism of clinically-utilized anti-cancer agents, disruption of microtubule networks may indicate a potential therapeutic mechanism of intense THz pulses. are the gene expression of and microtubule-associated genes in and the predicted dysregulation of cytoskeleton-related processes at the of biological

. Microtubule structure and fluorescence microscopy image. (a) Bound αβ-tubulin dimers comprise the hydrogen-bonded protofilaments [8]. 13 protofilaments laterally bind to form a hollow cylindrical polymer with a 25 nm diameter. Green outlines indicate a tubulin dimer that comprises a protofilament, and a single protofilament that comprises the polymerized MT. (b) Fluorescence image of rhodamine-labelled polymerized microtubules in solution, with an average length of 10 µm. The inset shows a zoomed window of polymerized MTs in a 20 µm square window. been observed to induce significant effects at multiple levels of biological organization that are not explained by the estimated heating [15][16][17][18]. While the specific nature of the fundamental interactions is not yet understood, these phenotypic effects may be attributable to interaction with sub-cellular macromolecular structures such as DNA [19][20][21][22] or cellular membranes [23][24][25]. However, experimental investigation of THz exposure effects to cytoskeletal protein complexes is still at a preliminary stage. Studies by Yamazaki et al. show that THz exposures of relatively low intensity (∼0.6 W/cm 2 ) can stimulate polymerization of actin, while highly intense THz pulses (∼3 GW/cm 2 ) resulted in actin disassembly [26,27]. To date, no studies have directly investigated the structural changes to MTs irradiated by intense THz pulses, although there have been indications in the literature that this may be an important interaction mechanism [28][29][30][31]. Hintzche et al. observed disturbance to the spindle apparatus (a mitotic structure comprised of microtubules) in hybrid animal cells [28], and Amicis et al. observed significant increase of micronuclei in THz-exposed human epithelial cells, which arise due to abnormal mitotic activity forming multiple nuclear envelopes in non-viable daughter cells [29].
Simulation studies investigating MT structure in nanosecond-duration pulsed electric fields have been reported, and show significant structural dysregulation for field strengths in the range of 50-750 kV/cm [32][33][34]. These studies predict reduction of MT stability proportional to field strength through conformational changes to key loops involved in lateral protofilament contact ( Fig. 1(a)), or alteration of local electrostatic properties at the GTP binding site. These results theoretically motivate the hypothesis that intense picosecond-duration pulses may also induce similar field-induced structural effects; however, more work simulating tubulin/MT response specifically to intense picosecond-duration fields is necessary. Simulation studies of THz fields for other macromolecular structures, such as membranes, offer insight into the expected interaction with intracellular environments [23]. For long-duration pulses (>∼10 ns), membrane interactions are predominantly due to ionic redistribution, as the field variations allow sufficient time for these motions to occur [35,36]. In contrast, for fast, short-duration pulses (<∼1 ns), the interaction is predominantly dielectric, arising from faster dipole reorientation dynamics. Importantly, this suggests that picosecond-duration pulses with sufficient field strength may interact with interior cellular structures, since the electronic screening effect in the membrane does not occur [35].
In this paper, we present the first observations of THz-induced disassembly of polymerized and chemically stabilized microtubules. Within minutes of exposure to a train of broadband intense THz pulses, we observed intensity-dependent MT disassembly. Further, bandpass-filtered exposures suggest additional frequency-dependence. These results assist in identifying and characterizing a fundamental interaction mechanism of THz radiation with biological systems that may underlie observed phenotypic responses, which is currently considered an open question.

Growth of fluorescently labelled microtubules
Rhodamine-labelled fluorescent tubulin (excitation/emission at 543/590 nm) was obtained from Cytoskeleton (Denver, USA), and MTs were polymerized according to the manufacturers' instructions [37]. One aliquot of rhodamine tubulin (Cytoskeleton, TL590M) was resuspended to 5 mg/mL (45.5 µM) in 4.0 µL of general tubulin buffer (Cytoskeleton, BST01) supplemented with 0.01 volumes of 100 mM GTP (Cytoskeleton, BST06) and 1 µL of MT cushion buffer (Cytoskeleton, BST05) to facilitate polymerization. This solution was mixed with 70 µL of an equivalent solution using unlabeled tubulin (Cytoskeleton, T240) for a final labelling ratio of 1:15. Aliquots were placed in an incubator at 37 • C for 20 minutes to allow the tubulin to polymerize to MTs of appropriate length (average ∼10 µm). Upon polymerization, the MTs were removed from the incubator and stabilized with 100 µL of taxol/MT buffer solution (20 µM taxol in general tubulin buffer). These were stored at room temperature as recommended, and the MTs were verified to be stable in this state for several hours. Further dilutions of this stock solution with taxol/MT buffer were performed as experimentally required. A fluorescence image of rhodamine-labelled polymerized and stabilized MTs is shown in Fig. 1(b).

Generation and detection of intense terahertz pulses
Intense THz pulses were generated by optical rectification of tilted-pulse-front laser pulses in lithium niobate (LiNbO 3 ), as schematically shown in Fig. 2(a). An oscillator/amplifier (Coherent Micra/Legend) generated a 1 kHz train of 800 nm, 50 fs, 3.6 mJ infrared laser pulses. An 1800 mm −1 reflective diffraction grating (RDG) established a pulse-front-tilt in the pump laser to satisfy velocity-matching conditions in the LiNbO 3 generation crystal with a 63 • cut output face [38]. Two cylindrical focusing lenses (f L1 =100 mm, f L2 =60 mm) in a 4f configuration imaged both the laser pulse front and grating surface onto the crystal output for optimal THz emission, as described in [39]. The intense THz pulse source for tilted-pulse-front optical rectification in lithium niobate (LN, LiNbO 3 ), using an 1800mm −1 reflective diffraction grating (RDG) and a pair of 4f-imaging lenses (L1 and L2), as described in [39]. Cross-absent bandpass filters (BPF) are used to isolate individual frequency bands. The THz beam is focused to either the sample location (beam propagating in the + z direction out of the page), or to the EO sampling system (beam in the -x direction) with a rotating gold off-axis parabolic mirror (ROAPM). A fraction of the pump pulse energy from a beamsplitter (BS) is attenuated (AT) and propagated colinearly with the THz beam for EO sampling in gallium phosphide (GaP). (b) The ROAPM is set at 0 • (+z, upwards) for through-substrate exposure of MTs in solution. The fluorescence excitation line is focused through the hole in the mirror and propagates to the sample colinearly with the focusing THz beam. The sample fluorescence emission is collected by long-working distance objectives, passed through a 578±16 nm bandpass filter, and analyzed in real-time with a CCD camera. For EO sampling, the ROAPM is set to 90 • .
The collimated THz emission from the LN crystal was filtered with black polyethylene to remove pump leakage. The beam was magnified and re-collimated by a pair of gold off-axis parabolic mirrors (OAPM) with focal lengths of f OAPM1 =15 mm and f OAPM2 =101.6 mm for a total magnification of 6.8. The expanded, collimated beam is directed towards a 76.2 mm focal length OAPM mounted to a programmable rotation stage (Thorlabs, PRM1Z8), as shown in Fig. 2(b). This provided 360 • control of the THz focus in a 76.2 mm radius annulus to focus to either the sample location (0 • ), or to the electro-optic (EO) sampling module for waveform detection (90 • ).
Temporal waveforms were detected by free-space EO sampling in a 200 µm (110) gallium phosphide (GaP) crystal mounted to a 2 mm (100) GaP substrate. Peak electric fields, E THz , were calculated directly from the normalized EO signal, ∆I/I 0 , from balanced photodetectors (BPD) according to the equation, where λ 0 =800 nm is the vacuum wavelength, n 0 =3.18 is the nominal GaP index, t GaP =0.46 is the GaP transmission coefficient, and r 41 =0.88 pm/V is the GaP EO coefficient [39]. Here, ∆I is the single-channel modulation of the BPD, and I 0 is the reference signal with no THz field present. The pulse energy and transverse intensity distribution were measured with a pyroelectric detector (Spectrum Detector, SPJ-D-8) and camera (Electrophysics, PV320), respectively.
Focused spot sizes were calculated as the 1/e 2 width of the Gaussian fits to 1D line profiles, and pulse durations were determined as the 1/e width of the waveform Hilbert transform. For investigations of frequency dependence, THz bandpass filters (Thorlabs, FB19M590/FB19M200, 19.6 mm inner diameter) were used as additional filtration to isolate the 0.5 THz and 1.5 THz bands. The broadband and bandpass waveforms, power spectra, and spatial intensity distributions at the focus are shown in Fig. 3. Typical pulse parameters are summarized in Table 1.

Alignment procedure for real-time analysis
To ensure the THz focus, sample location, and microscope imaging plane were coincident in space, the alignment procedure in Fig. 4 was employed. First, the THz pulse energy was maximized through a 1 mm pinhole alignment aperture placed at the sample location with a pyroelectric detector (ScienceTech, 6925-01) centered above. Once this signal was maximized, the detector was removed, and the microscope field-of-view (FOV) was focused to the pinhole to ensure the imaging plane and THz focus were longitudinally aligned and centered. Next, a sham sample was loaded into the holder and longitudinally translated until the sample was in focus. This achieved coincidence of the sample and image planes with the THz beam focus for the real-time fluorescence analysis of MTs.  1) The dish is levelled, and the THz beam focus is located by translating the holder and maximizing pulse energy through a 1 mm pinhole aperture centred over the beam input window. (2) The microscope is centred (x-y) and longitudinally (z) aligned to the THz focus by bringing the alignment aperture used for THz energy measurements into focus in the brightfield FOV.
(3) The sample plane for a given substrate is aligned to the THz focus by longitudinally (z) translating the sample holder until a test sample is in focus.

Microtubule exposures
MT solutions were loaded onto double-well microscope slides of optical plastic that were verified to be transmissive (T=0.92) to THz wavelengths (Ibidi µ-slide, Cat. 80281), with sufficient volume to fully cover the focused THz spot once compressed with a coverslip (1-3 µL). The MTs were left for several minutes to 1 hour to settle to the substrate and ensure stability. For exposure, the THz beam was focused through a 5 mm diameter hole in the sample holder with the rotating OAPM to one of the wells, and the second well provided matched, unexposed control solutions. Similar to the recommended MT stock storage conditions, exposures were carried out at room temperature.

Real-time fluorescence imaging and image analysis
THz-exposed rhodamine-MTs were imaged in real-time with a fluorescence microscope aligned to the beam focus. The fluorescence excitation laser (532 nm CW diode-pumped solid-state laser) was focused from below through a 1 mm hole in the rotating OAPM, and uniformly excites the sample growth area. Short pass filters eliminated any residual 805 nm and 1064 nm emissions, and a neutral density filter wheel tuned the incident power. A programmable shutter connected to the imaging camera via USB blocked the laser and allowed fluorescence excitation to be controlled within the imaging software. Fluorescence emission (∼590 nm) from the sample plane was collected by long-working distance 5x/20x objectives (Mitutoyo, 378-802-6/378-804-3) and a 2x tube lens (Edmund Optics, MT-2, 56-863), passed through a 578±16 nm bandpass filter, and collected by a 2.8 MP monochrome CCD microscopy camera (Lumenera, Infinity 3-3URFM).
For each timepoint acquired during THz exposure, a series of 15 images were averaged with an exposure time of no greater than 300 ms each. The excitation laser therefore exposed the samples for no more than ∼5 s for each image set. Excitation control images were acquired, and a photobleach calibration curve was experimentally determined. It was verified that the timescale for the fluorescent label to experience significant levels of photobleaching was much greater than the excitation exposure times utilized in this study. For imaging, auto-exposure control (AEC) was enabled when investigating general structural effects. In cases where the pixel intensity values were utilized for quantitative calculations (Section 3.2), AEC was disabled. Sample fluorescence images for low and high tubulin concentrations are shown in Fig. 5. Low initial tubulin concentrations result in individually resolvable MTs, while large tubulin concentrations form large aggregate structures. Image processing and analysis was performed using the open-source software, ImageJ (https://imagej.net). Fluorescence images were background-subtracted with a rolling-ball algorithm. For the frequency-dependent studies, the structural change to MTs were quantified by the area fraction of the image space occupied by MTs. MT area fractions were determined by first creating a binary image mask for each timepoint with a common pixel intensity threshold.
As MTs disassemble, intensity per unit area reduces, which reduces the total fraction of pixels with rhodamine fluorescence signal intensity above the chosen threshold for a given imaging region. The MT area fractions were determined with the ImageJ plugin "Analyze Particles" applied to the masked images, which algorithmically contours the masked MTs and determines the enclosed areas. Each set of MT area fractions were normalized to the initial timepoint.

Results
3.1. Intensity-dependence of broadband terahertz pulses on microtubule structure Figure 6 shows three sets of time-series images from broadband THz exposures for varying exposure conditions. Differing combinations of imaging magnification (M) and tubulin concentration (C tub ) were utilized to analyze either single or aggregate MT structures in either uniform or non-uniform intensity distributions.
High-M and low-C tub experiments (40x, 0.25 mg/mL) provide detailed structural resolution to individual MTs. As shown in Fig. 6(a), disassembly of a single MT is shown to occur within 11 minutes with 0.8 µJ pulses and a 230 kV/cm peak field strength. A multimedia file of the single MT breaking in Fig. 6(a) is included (see Visualization 1).
Low-M and High-C tub (10x, 5 mg/mL) experiments with similar exposure parameters show the large-scale effects to MT aggregates ( Fig. 6(b)). The large FOV allows analysis of intensitydependence within a single exposure: By analyzing structural effects to MTs for varying distances from the center of the THz focus, differential effects to MT structure for varying energy densities are investigated. From within the total imaging FOV of 870×655 µm 2 , three 100×100 µm 2 FOVs were selected for analysis, corresponding to different regions of the THz intensity distribution (shown at right in Fig. 6(b)). In the central region, the energy density reaches 80 µJ/cm 2 , and the largest qualitative change to MT structure is observed. Near the beam edge ∼0.5 mm from center, the energy density falls to ∼30 µJ/cm 2 , and no significant structural change is observed.
High-M and High-C tub (40x, 5 mg/mL) experiments show structural effects to MT aggregates in a nearly uniform THz intensity profile at higher THz pulse energy (1.2 µJ) and field strength (400 kV/cm and 409 kV/cm), depicted in Fig. 6(c). In both time-series, MT disassembly and aggregate destruction occurs significantly faster than the previous exposure cases, within ∼5 minutes, further corroborating the evidence for intensity-dependence that was indicated in the differential effects seen in the low-magnification images.

Exposures with terahertz bandpass filters
From the broadband exposure results in Section 3.1, it is not possible to conclude that the observed effects are strictly only intensity-dependent, since wavelength-dependent focusing induces spatial variation in frequency for the broadband pulse, in addition to the power-density variation. The central intense region of the beam contains a larger fraction of high frequency energy relative to the low-intensity beam regions, and so the differential effect to large MT aggregates for varying location in the THz focus may additionally be influenced by the frequency distribution associated with the region of space under consideration.
To isolate the potential effects for varying spectral content, THz bandpass filters are used to transmit narrow 0.5 THz and 1.5 THz bands, and compared to broadband exposure, adjusted for intensity differences. For these analyses, THz-induced MT disassembly is quantified using the MT area fraction. As shown in Fig. 7(a), contours determined with the analysis algorithm (Section 2.5) define this area fraction, which decreases over time as the fluorescence signal from disassembled MTs falls below the mask threshold. The results in Fig. 7(b) show that the decrease in MT area generally follows an exponential decay. For comparison, an equivalent analysis with images of unexposed MTs is included. The curves are labelled with characteristic times from fitting these data to an exponential function to establish quantitative timescales of THz-induced effect for varying spectral content. The top plot of Fig. 7(b) shows the total pixel intensity of the images for each timepoint. The maximum loss of total fluorescence signal was ∼5%, and so the decrease in MT area fraction of THz-exposed samples is not due to photobleaching of the rhodamine label, but rather to MT disassembly, as illustrated in the second column of Fig. 7(a). By converting quantitative fluorescence images (column 1) to a binary image with a common threshold (column 2), the disassembly of MTs over time may be quantified by the change of area fraction with rhodamine signal above a common intensity threshold. The area fractions are determined by algorithmic contours with the ImageJ plugin, "Analyze Particles". The reduction of area fractions follow an exponential decay curve. (b) Fractional MT area calculated using the procedure in (a) for varying THz bands, with dashed curves representing exponential fits, and τ is the associated characteristic time. Each dataset was normalized to the initial relative MT area. The exponential fit qualities (R 2 ) are 0.99, 0.96, and 0.92 for the broadband, 0.5 THz, and 1.5 THz fits, respectively. An equivalent analysis on unexposed MTs is included for reference. Top: The total pixel intensity of the raw images. The total rhodamine signal of all images does not degrade significantly (<5%) over the exposure duration, indicating the MT area fraction decay is not due to photobleaching. (c) The MT area fraction vs. total dose (J/cm 2 ), which corrects for differences in pulse energy and focused spot area (see Table 1). D X is the characteristic dose for the corresponding curve (e −D/D X ) . The characteristic total dose for the low-frequency 0.5 THz band (1.4 J/cm 2 ) is significantly lower than both the broadband and high-frequency 1.5 THz band (13 J/cm 2 and 23 J/cm 2 , respectively), indicating frequency-dependence of THz-induced MT disassembly, with greater disassembly induced by low-frequency THz energy (∼0.5 THz).
While there is a large range of characteristic times for varying spectral content from the fits in Fig. 7(b), this could be due to differences in pulse energy and focused spot area between the separate THz pulses (Fig. 3, Table 1). To correct for these differences, the MT area fraction is plotted as a function of the total dose, D = I avg t, as shown in Fig. 7(c). The total dose is reported in units of J/cm 2 , and is distinct from the per-pulse energy density (µJ/cm 2 ) in Fig. 6. The difference in characteristic doses indicate that THz-induced disassembly of MTs is significantly influenced by the pulse frequency content, as discussed in the following section.

Discussion
Exposures at higher magnification (40x) and low tubulin concentration (0.25 mg/mL) provide an opportunity to observe detailed effects resolvable to single MTs (Fig. 6(a)). Conversely, analysis of low-magnification (10x), high-concentration (5 mg/mL) MT samples are used to investigate intensity-dependence in a single exposure by analyzing separate regions of the imaging FOV for varying locations in the Gaussian intensity distribution (Fig. 6(b)). In the former case, disassembly of a single MT polymer is observed to occur within 11 minutes of THz exposure. In the latter case, significant variation in THz-induced change to MT structure is observed for varying locations within the THz focus over the exposure duration. The strongest effect is observed in the central region with the highest intensity (top row of Fig. 6(b)), in which the blurring of large aggregate structures is attributed to single MTs disassembling, as in Fig. 6(a). As the analysis window moves to the lower-intensity outer area of the beam, the magnitude of the effect diminishes.
The above cases were exposed with a moderately attenuated THz beam (0.8 µJ/pulse, 230 kV/cm), and the disassembly to MTs occurred in 11 minutes. Figure 6(c) shows high-concentration exposures with higher pulse energy (1.2 µJ/pulse) and field strength (400 kV/cm and 409 kV/cm), and at high magnification such that the THz intensity profile is roughly uniform across the imaging FOV. Here, initiation of MT structure changes is observed within 4-6 minutes, further supporting the intensity-dependence observed above.

Frequency-dependence
Exposures at high magnification (40x) and high tubulin concentration (5 mg/mL) with THz bandpass filters investigated coarse frequency-dependence within an approximately uniform intensity distribution across the imaging FOV. The characteristic times for MT disassembly extracted from the exponential fits in Fig. 7(b) are 8±2 min (R 2 =0.9626) and 22±6 min (R 2 =0.9223) for the 0.5 THz and 1.5 THz bands, respectively. While these are significantly longer than the broadband reference case (τ=3.1±0.6 minutes, R 2 =0.9854), they also correspond to significantly lower pulse energy, intensity, and peak electric field, as summarized in Table 1. Figure 7(c) shows the MT area fraction curves in terms of dose delivered (J/cm 2 ), which corrects for differences in pulse energy and area, but not peak electric field. If effects depended only on intensity, the corrected curves are expected to overlap. However, the characteristic dose for the 0.5 THz band is 1.4±0.3 J/cm 2 , relative to 23±6 J/cm 2 for the 1.5 THz band. Additionally, the field strength of the 0.5 THz band (37 kV/cm) is significantly lower than the 1.5 THz field strength (100 kV/cm). Since the 0.5 THz disassembly occurs faster not only at a lower dose, but also a lower peak field strength, these data indicate that MT disassembly is significantly influenced by frequency content. The broadband curve corresponds to a characteristic dose of 13±3 J/cm 2 , intermediate to the curves for the low and high frequency bands. This suggests that the observed broadband MT disassembly is predominantly due to low-frequency energy, while much of the high-frequency content is not utilized for MT disassembly, but is likely instead absorbed by the media or thermalized.
Other effects may potentially explain the observed frequency dependence as well, such as differential absorption in the aqueous media environment prior to MT interaction. In water, the absorption coefficients at 0.5 THz and 1.5 THz are approximately 82.6 cm −1 and 153.6 cm −1 , respectively, and so more of the high-frequency band will be absorbed in aqueous media rather than the MTs [40]. Additionally, in MTs the calculated absorbed energy at 1.5 THz is also roughly twice the energy absorbed at 0.5 THz, and so these competing effects obscure one another experimentally [5]. Future studies investigating concentration dependence (i.e., varying probabilities of THz absorption in MT structures relative to surrounding media) or variation of aqueous media layer thickness (i.e., varying interaction distance in the MT solution) will assist in contextualizing these apparent frequency-dependent effects to better understand the precise nature of this interaction.

Considerations of thermal or shockwave interaction mechanisms
An important consideration is the potential for other interaction mechanisms, such as thermal effects or field-induced shockwaves, to influence the biological response in THz exposure studies [15,27]. In our experiments, the duty cycle is limited (1 kHz train of picosecondduration pulses) to ensure negligible heating in the biological media. An estimate of the maximum per-pulse temperature increase from the broadband THz pulse is roughly 5 mK, using ∆Q = mc∆T = ρVc∆T, where ∆T is the temperature change due to a pulse energy ∆Q, and assuming similar properties to water (density ρ=1 g/cm 3 , specific heat capacity c=4.2 J/g/K). The volume V = Az is calculated from the THz spot area and penetration depth. The average steady-state heating due to the 1 kHz pulse train is measured with a thermal imager (Reed Instruments, R2100) to be less than 1 • C in water, and is consistent with measurements of other similar THz exposure systems [41][42][43]. These are also consistent with the THz-water heating model from Kristensen et al., which predicts a maximum temperature change at the beam center of 2 • C for our broadband beam parameters [44].
In addition to the measurements above, experiments with taxol-stabilized MTs inherently control for thermal effects, as they are verified to be stable for several hours at room temperature, and higher temperatures assist MT polymerization (indeed, it is the mechanism exploited to polymerize MTs from tubulin by placing them into a 37 • C incubator, as described in Section 2.1). Therefore, disassembly by THz represents the opposite effect to that expected from increasing sample temperature. This implies that the dominant interaction mechanism is non-thermal and may instead be explained by coupling to natural oscillatory dynamics of MT structures.
Using intense THz pulses generated by a free electron laser, Yamazaki et al. induced dramatic disassembly of actin polymers within 30 minutes that was attributed to an acoustic shockwave formed in the aqueous medium that penetrates significantly deeper than EM THz energy [27]. While this represents an exciting new interaction mechanism to explore that may play a role in THz-induced biological effects, we cannot claim with confidence that similar shockwave phenomena are occurring in the present study. Tsubouchi et al. show that the acoustic shockwave amplitude is expected to scale linearly with the product of absorption coefficient and fluence, αF, with significant shockwave amplitude requiring αF on the order of ∼0.1-1 J/cm 3 [45]. The relatively lower values in our work (αF < 0.02 J/cm 3 ) may not be sufficient to produce acoustic waves with significant amplitude, although this should be explored as a potential mechanism in future investigations. Nevertheless, we observed MT disassembly within minutes at similar fluence, but with lower acoustic generation efficiency.

Relation to tissue-level biological effects
Recently, we have reported on differential gene expression induced by intense THz pulses in human skin, and identified biological processes predicted to be dysregulated [46]. In these experiments, 3D human skin tissue models were exposed using similar beam parameters (2.4 µJ/pulse, E peak = 240 kV/cm, I pulse = 74 MW/cm 2 ) as those outlined in Section 2.2. Here, we emphasize additional aspects of these data in the context of the MT exposures presented in this paper. Specifically, these are the THz-induced gene expression of tubulin and microtubuleassociated genes in human skin, and the predicted dysregulation of cytoskeleton-related processes at the tissue-level of biological organization. Figure 8(a) shows the differential gene expression induced by intense THz pulses of all genes in the tubulin superfamily and genes that encode for microtubule-associated proteins (MAPs). Intense THz pulses largely downregulate expression of tubulin and other genes related to MT structure and function. There is a particularly large suppression of genes in the α and β subfamilies (TUBA and TUBB), which encode for the principal structural components of MTs, as described in Section 1. Gamma and delta tubulin, which play critical roles in MT nucleation, were unaffected. In the context of the significant disassembly of MTs induced by THz pulses presented in this study, this differential gene expression may represent the cell's genomic response to a disassembled and disrupted cytoskeletal network.
Downregulation of TUBA and TUBB following chemical disassembly of MTs in human cells has been established for decades [47][48][49][50], further corroborating the significance of THz-induced MT disassembly as a key mechanism driving the tissue-level biological response in skin in Fig. 8(a). However, more work needs to be done to elucidate the precise nature of this potential THz interaction mechanism. Among the other microtubule-regulating genes affected by THz are the microtubule-associated protein (MAP) family that regulate MT growth, and PLK4, a serine/threonine kinase that regulates mitotic centriole dynamics [51]. With the exception of TTL, a cytosolic enzyme responsible for post-translational modifications of alpha tubulin, all significant THz-induced effects on MT-associated gene expression were inhibitory.
Some of these THz-affected genes have also been investigated for therapeutic application in gene therapy: Knockdown of class II, III, or IV β-tubulin family genes (TUBB2A, TUBB2C, TUBB4Q, and TUBB3 in Fig. 8(a)) has been shown to enhance the effectiveness of tubulinbinding pharmacological agents via suppression of MT dynamics and sensitization to apoptosis induction [52,53]. These data motivate the hypothesis that intense THz pulses may enhance the sensitivity of diseased cells to similar drugs via interaction with MTs and potentially other cytoskeletal structures.
Gene Ontology (GO) analysis was performed on the measured expression dataset to identify statistical over-representations of biological processes, cellular components, and molecular functions among the set of significantly differentially expressed genes, relative to expected global background rates [54]. 58 total GO terms were identified as over-represented by the THz-induced expression profile (p<0.01). These were largely related to epidermal processes (as the sample under study was skin), however, 8 GO terms were associated with general or specific cytoskeletal processes, components, and functions. The GO terms, odds-ratios (OR), and p-values are shown in Fig. 8(b), and include over-representation of THz-affected genes in several important cytoskeleton-regulated processes. Thus, intense THz pulses may dysregulate these cytoskeletal-related functions leading to phenotypic changes observed at cellular or tissue levels. In the future, studies explicitly investigating THz-induced changes to MTs or other cytoskeleton protein complexes correlated to the expression dynamics of genes identified in Fig. 8(a) will assist in clarifying the precise nature of the cellular response to THz-MT interactions.
The present study shows that intense THz pulses are capable of dissociating cellular concentrations of polymerized, taxol-stabilized MTs within minutes of application, and are consistent with differential gene regulations measured in cellular systems. The onset of MT disassembly is dramatic and clearly triggered by THz power input, as seen in repeated, controlled experiments. Nonetheless, several exposures have also been performed in which no effect to MT structure was observed, even in similar exposure conditions that had previously induced disassembly. Since MTs are verified to be stable at room temperature for hours, spontaneous disassembly in controlled exposure studies is likely not occurring. The most likely explanation is that our current THz sources, while sufficient to induce MT disassembly under optimal conditions, are close to the levels necessary for detectable MT disassembly, as suggested by the diminished effect near the outer regions of the beam, but still within the defined THz focus (Fig. 6(b)). Small variations in the MT environment may significantly affect the THz energy and field such that the pulse seen by the sample is rendered insufficient for significant MT disassembly. These data suggest that THz-induced disassembly requires energy density in the range of 30-80 µJ/cm 2 , a peak-field strength > 37 kV/cm, total dose on the order of ∼1-10 J/cm 2 , and a further dependence on spectral content of the THz pulse, with faster disassembly for low-frequency THz bands (∼0.5 THz). More research investigating the THz parameter space for differential structural change is necessary.
The inhibition of MT dynamics is a standard mechanism of action for several types of pharmacological agents used in chemotherapy. Paclitaxel (also known as taxol, the MT stabilizing agent used in these experiments) is an early chemotherapy drug that inhibits MT dynamics by stabilizing polymerized filaments such that necessary cycles of polymerization/depolymerization ("treadmilling") are suppressed [3,55]. Similarly, colchicine is an experimental anti-cancer (and anti-inflammatory) drug that inhibits tubulin dimers from polymerizing into MT filaments [56]. Several other types of tubulin-binding pharmacological agents such as vinblastine, demecolcine, or nocodazole, similarly inhibit MT dynamics and reduce polymer mass [57].

Conclusions
In this paper, microtubule disassembly by intense THz pulses was investigated. Significant disassembly of MTs is observed within several minutes of THz power input, and these effects are not explained by heating or shockwave formation. Significant intensity dependence was observed in the rate of MT disassembly. Moreover, exposures using THz bandpass filters suggest that the frequency content has a significant influence on the detected structural changes. However, more work is required to isolate confounding variables in exposure experiments, such as the reduction in THz energy transmitted through aqueous media, relative to the increased absorption in MTs at higher frequencies. Our results were analyzed in the context of our previous THz exposures performed in skin tissue, and suggest that MT disassembly by intense THz pulses may be a key interaction mechanism driving the biological response observed in higher-level systems like multicellular tissue.
As MTs are essential mitotic structures and popular targets of anti-cancer therapies, the THzinduced effects to MT structure reported in this paper suggest a potential therapeutic mechanism of intense THz pulses, possibly in combination with existing pharmacological interventions. Future research investigating these effects in higher-level cellular and tissue systems is necessary to establish the prospective clinical feasibility.