Actively Driven Fluctuations in a Fibrin Network

Understanding force propagation through the fibrous extracellular matrix can elucidate how cells interact mechanically with their surrounding tissue. Presumably, due to elastic nonlinearities of the constituent filaments and their random connection topology, force propagation in fiber networks is quite complex, and the basic problem of force propagation in structurally heterogeneous networks remains unsolved. We report on a new technique to detect displacements through such networks in response to a localized force, using a fibrin hydrogel as an example. By studying the displacements of fibers surrounding a two-micron bead that is driven sinusoidally by optical tweezers, we develop maps of displacements in the network. Fiber movement is measured by fluorescence intensity fluctuations recorded by a laser scanning confocal microscope. We find that the Fourier magnitude of these intensity fluctuations at the drive frequency identifies fibers that are mechanically coupled to the driven bead. By examining the phase relation between the drive and the displacements, we show that the fiber displacements are, indeed, due to elastic couplings within the network. Both the Fourier magnitude and phase depend on the direction of the drive force, such that displacements typically propagate farther, but not exclusively, along the drive direction. This technique may be used to characterize the local mechanical response in 3-D tissue cultures, and to address fundamental questions about force propagation within fiber networks.


INTRODUCTION
The extracellular matrix (ECM) interacts with and provides a variety of mechanical and biochemical cues to cells [1,2]. Experimental systems using naturally derived ECMs more closely mimic in vivo cell-ECM interactions than, for example, plastic [3]. Common naturally derived ECMs such as collagen and fibrin hydrogels exhibit important characteristics such as micron-scale pores and cell adhesion sites [4]. The fibrous nature of these ECMs plays potentially critical roles in cell-cell communication through the ECM, especially in 3-D culture systems [5,6]. For example, cells respond to forces an order-of-magnitude farther away in naturally derived fibrous materials than materials with sub-micron pores, often treated as a continuum at the cellular scale [7][8][9][10]. Therefore, understanding force propagation through naturally derived fibrous materials at the cellular scale is of great interest. Fibrin, a major component of hemostasis, is well suited for studying force propagation through fiber networks. Fibrin forms the provisional ECM of a blood clot, is used as tissue glue, and is a substrate for the study of wound healing and in tissue engineering [11][12][13][14]. Fibrin hydrogels can be formed by activating soluble fibrinogen with thrombin, which forms a randomly arranged 3-D network of fibrin fibers.
At the fiber scale, experimental observations using optical tweezers and atomic force microscopy (AFM) as well as computational simulations show that individual fibers act as strain stiffening elastic beams [15][16][17][18][19][20]. While force transmission at the scale of single fibers is well understood, force transmission in assembled networks has not been fully elucidated [21]. At scales much larger than the average pore size, the material can be treated as a continuum, but at the scale of pores, and even cells, it is necessary to consider that forces cannot elastically transmit through the liquid region and are confined to the fibers [9,22]. At this mesoscale, the existence of a preferred pathway for force transmission or a 'force transmission highway' has been observed in simulations and inferred in experimental observations [7,9,21,[23][24][25]. Simulations that take into account the fibrous structure of the ECM show that when the matrix is extended, subsets of fibers within the network carry higher tension than would be expected from a continuum model of the material, while other fibers carry less tension or are buckled under compression [22,23]. How and where these high tension pathways form and change may play a role in cell orientation [26], cell migration or invasion [27,28], and capillary sprouting [29]. A method for directly assessing which fibers in a network carry tension can provide insight into force transmission at this mesoscale.
Here, we describe an approach to detect which fibers within a 3-D fibrin hydrogel are involved in force transmission. Optical tweezers are used to apply a low frequency oscillatory force to a bead embedded in a fibrin network. We measure the response of the fibers within the network using confocal microscopy by analyzing the fluctuation of the fluorescence intensity. Both confocal microscopy and optical tweezers are able to access regions deep within a hydrogel without invasive disruption of the structure, unlike AFM. Future studies will aim to extend this method for identifying force transmission pathways around cells and to study how naturally derived ECMs participate in cell-cell and cell-ECM communication.

Bead Oscillation
A 1064 nm continuous wave laser (IPG, YLR-5-1064-LP) was used to optically trap a 2 m diameter silica bead embedded approximately 35 μm above the glass in a fibrin network with a 60×1.45 NA oil immersion microscope objective (Olympus). The focused laser beam was steered as a sine wave in the image plane using galvanometer mirrors (Thorlabs, GVS102). Bead displacement was detected by a co-aligned 785 nm (Thorlabs, LP785-SF100) laser-diode beam. The beam was back scattered from the bead and directed onto a quadrant photodetector (Newport, 2,901 and 2,903). Trap stiffness was determined by trapping a bead without oscillation of the trap and measuring the corner frequency of the displacement power spectrum in water [30]. Laser power was adjusted to achieve a trap stiffness of 34 pN/m for all experiments. Schematic of the optical tweezer system can be found in Supplementary Figure S1. Optical tweezers were used to oscillate the bead in the x (horizontal in the image) or y (vertical in the image) directions at the desired frequency (driving frequency, F d ). The position of the optical tweezers oscillated as a sine wave with 1 m peak-to-peak amplitude unless otherwise stated.

Fluorescence Microscopy
Fluorescence microscopy images were recorded using an Olympus Fluoview 1,200 laser scanning confocal microscope. The Atto 488 NHS-ester fluorophore was excited using the 488 nm laser line operating at 0.7% of maximum power as indicated in the Olympus Fluoview software. Images of fibrin were acquired using the rapid-scanning imaging mode to achieved maximum frame rates of 15.4 frames/second. This imaging mode limits image resolution to 256 by 256 pixels to cover a field of view of 21.1 μm by 21.1 μm unless otherwise stated. 400 consecutive images were recorded over 26 s. Z-stacks of at least 10 μm thickness (5 μm above and below the bead) were acquired prior to bead oscillation in order to identify the local fiber network around each bead.

Image Analysis
Images in the Olympus OIF format were converted to TIFF files using ImageJ [31]. Pixel-by-pixel 1-D Fourier transforms were performed in MATLAB (MathWorks). Figures were generated with MATLAB. Binary image masks were generated by setting a threshold of mean intensity (in time) of each pixel that is 20% greater than mean intensity (in space and time) of all pixels. Pseudo-color 3-D image stacks were generated using the 'temporal-color code' function in ImageJ and all figures were arranged in Microsoft PowerPoint. Our custom MATLAB code is made available in Supplementary Materials along with all TIFF images. One must first execute pixelflucatuation.m for image input and fast Fourier transform computation. Next, pixelflucatuationfigures.m should be executed. It operates on the output from pixelflucatuation.m to generate all figures. The code as provided assumes all files are placed in the same folder.

Quantification of Local Phase Variance
The Orientational Order Parameter (OOP), also known as the nematic order parameter [32,33], was used to determine local organization in phase (θ) images. For OOP, the variance in the orientation of n unit pseudovectors ( r → , which has components r x cosθ and r y sinθ) can be quantified [34][35][36]. Pseudovectors are vector-like objects, which are invariant under inversion ( r → − r → ), and are also called "axial vectors" [37]. The OOP ranges from 0 to 1, and will be 0 when the data is uniformly distributed (perfectly isotropic) and will increase as the data becomes more anisotropic toward a value of 1 for perfectly aligned data.
In order to calculate the OOP, the tensor, The orientational order tensor, T in Equation 2, is constructed by normalizing each tensor T i (Equation 1), and averaging the normalized tensors. The OOP is the maximum eigenvalue of the orientational order tensor T, shown in Equation 2 T 2 r i,x r i,x r i,x r i,y r i,x r i,y r i,y r i,y − 1 0 0 1 (2) In order to quantify local phase organization, the OOP was calculated at each pixel in a 7 × 7 pixel neighborhood. The MATLAB functions calculate_OOP.m and hoodOOP.m, found in the Supplementary Materials, were used for the OOP calculations.

RESULTS
Fluorescence confocal microscopy videos were recorded while optical tweezers applied an oscillatory force on a 2 μm bead entangled in a fibrous fibrin network. Supplementary Video S1 shows a bead with no applied oscillation. Supplementary Video S2 shows the same bead oscillated in the horizontal x-direction. The bead and some of the surrounding fibers move noticeably in response to the oscillation of the optical tweezers. For many fibers, motion is difficult to track precisely from frame to frame due to noise. However, a subset of pixels show increased intensity fluctuations at the optical tweezers' drive frequency (F d ) ( Figure 1). Intensity, denoted as I, is plotted for several representative pixels overlaying the bead, fibers, or liquid region ( Figures 1A,B). Figure 1A shows the mean intensity image, where each pixel is the mean fluorescence intensity over time, denoted by < I >. The spectral magnitudes of the 1-D Fourier transform of I in time, denoted by I , for the representative pixels are plotted in Figure 1C. I does not show prominent peaks without applied oscillations ( Figure 1C orange). However, with applied oscillations, some of the representative pixels show a peak at the driving frequency, denoted by F d (Figure 1C blue). Pixel i lies on the bead and a strong peak in I occurs at F d . A similar but less prominent peak is present for pixel ii, which lies on a fiber directly attached to the bead. Pixel v does not lie on a fiber and is instead in the liquid region or pore of the fiber network. We found that I in these regions are insensitive to the applied bead oscillation and pixels in these regions also have lower I as well as lower I at all frequencies. Pixels iii and iv lie on fibers farther away from the bead, which are not directly connected to the bead by the in-focus portion of the fiber network ( Figure 1A). Determining if these fibers are oscillating with the bead can be difficult based on the spectral plot of individual pixels as presented in Figure 1C, when I at F d is near to the mean I across all frequencies. We can improve identification of fibers which are mechanically coupled to the bead by taking into account the intensity fluctuations of all the pixels over each fiber. We can visualize these spatial relationships in I images (Figure 2, Supplementary Video S3 bottom). As an example, Supplementary Video S3 shows the I images series for the same bead as in Figure 1 with (bottom right) and without (bottom left) applied oscillation. Even without applied oscillations, fibers are visible in these I images (example: Figure 2A) as regions of increased I relative to the liquid regions that form the pores within the fibrin network. When an oscillation is applied, I images are similar to the nooscillation condition, except at F d (Supplementary Video S3 bottom, Figures 2A-D) where regions of elevated I are evident and lie on a subset of fiber pixels. To remove the portion of I which is independent of applied oscillations, we first define a background image 〈 I 〉 to be the mean value of I across all frequencies, excluding 0 and F d ( Figure 2E) and then < I > is subtracted from each I image (example: Figures 2F-H). For the remainder of this paper we redefine I as I − < I > to simplify notation, unless otherwise noted. The I image at F d ( Figure 2G) highlights oscillating pixels, but cannot identify the fibers that do not oscillate, since fibers not oscillating in Figure 2G are indistinguishable from the liquid region. To identify both oscillating and non-oscillating fibers, image masks can be generated from < I > images (example: Figure 1A) and applied to I images (example: Figures 2I-K). Reexamination of pixels iii and iv from Figure 1 shows that the fiber containing pixel iii is oscillating in response to the bead while pixel iv is not ( Figures 2G,K). In summary, an increase in pixel intensity fluctuations can be detected for subsets of pixels when oscillatory forces are applied to the fiber network. Such pixel intensity oscillations preserve the frequency of the applied force and are identifiable after mean-magnitude subtraction.
In addition to I , we investigated the corresponding phase of pixel fluctuations, ∠ I (Supplementary Video S3 top). In the absence of applied oscillations, ∠ I is random and uniformly distributed ( Figures 3A,B). When an oscillation is applied to the bead, the ∠ I at F d displays local organization along some of the fibers. These fibers appear as two joined parallel lines shifted in phase by π radians ( Figure 3D). The phase distribution of all pixels becomes non-uniform exhibiting modes separated by π ( Figures 3D,E). This π shift is an artifact of imaging due to the motion of the fibers relative to the fixed position of pixels, such that as a fiber moves rightwards, the pixels on the right side of the fiber get brighter. At the same time, pixels on the left side of the fiber get dimmer. The relationship reverses as the fiber moves leftwards. This distinct pattern allows for easy identification of the oscillating fibers by eye. We can distinguish the background liquid region and non-oscillating fibers which both have random ∠ I distributions by applying the 〈I〉 mask to exclude the liquid regions just like we did with the I image ( Figure 3G). Notably, the fibers which have increased ∠ I organization in Figure 3G match the fibers which show increased I in Figures 2K, 3C,F. The organization of the phase of the pixel intensity fluctuations at F d can be quantified using the Orientational Order Parameter (OOP), which treat the phase vectors as pseudovectors with components equal to cos(∠ I) and sin(∠ I) ( Figure 3H) ( [11,32,44]). A locally disordered system will have a local OOP of 0 while a locally ordered (i.e. aligned) system will have a local OOP of 1. In this work, we define "local" to be pixels in a seven by seven pixel neighborhood. In summary, fiber fluctuations resulting from forced bead oscillations are detectable by ∠ I and I , while local organization of ∠ I can be quantified using OOP. We have shown that oscillation of the optical tweezers with an amplitude of 1 m results in a response from the bead and a subset of the fiber network detectable by the ∠ I, I , and OOP. Next, we investigate how ∠ I, I , and OOP change with oscillation amplitude or direction (Figure 4). There is a subtle increase in the I as the amplitude of oscillation was increased from 1 μm to 1.2 μm ( Figure 4A). While it is difficult to directly identify differences in the ∠ I in response to increasing the amplitude of oscillation, the increased fiber recruitment and phase organization of already oscillating fibers is apparent when looking at the OOP image ( Figures 4B,C). Conversely, decreasing the amplitude of oscillation results in decreasing I and OOP (Figure 4). At an oscillation amplitude of 0.2 μm, the ∠ I, I , and OOP are indistinguishable from the no-oscillation control (Figure 4). Subtle differences appear in the ∠ I and OOP images with oscillations at or above 0.5 m, but not the I image ( Figure 4). This indicates that the ∠ I is more sensitive to small oscillations in the fiber network than the I . Now that we have shown pixel fluctuations are dependent on the optical tweezers oscillation amplitude, we next investigated if the subset of image pixels that are phase coupled (i.e. on fibers) changes with bead oscillation direction. Experiments are repeated with forced oscillations in the y-direction after oscillations in the x-direction, both with 1 μm oscillation of the optical tweezers. The in-focus fiber network structure can be seen in the 〈I〉 image for each bead ( Figure 5A). It is observed that some subset of the fiber network responds to the oscillation of the bead and this response can been seen in the I , ∠ I, and OOP images ( Figures  5B-G). The subset of engaged fibers is dependent on the direction of the bead oscillation ( Figure 5 arrows). We can compare the response of the fiber network to the direction of oscillation by looking at the differences in I images between applying oscillation in the x-direction or y-direction. It is apparent that some regions of the fiber network are more responsive to oscillation in the x-direction (blue) or oscillation in the y-direction (red) ( Figure 6A). This preference in response suggests fluctuations travel farther through the network along the axis of oscillation. This trend can be visualized by comparing I of pixels located along the oscillation direction of the bead to the pixels located perpendicular to the oscillation direction ( Figure 6B). Four regions of the I image from Figures 5B,C, each π/4 radians wide, are denoted as either parallel (cyan) or perpendicular (magenta) to the axis of oscillation ( Figure 6B). Figure 6C shows that I of the fibers decreases with distance from the bead, but to a lesser extent in the direction parallel to the oscillation. The force transmission observed in individual experiments is complex and not easily predicted from the 〈I〉 image. For example, when bead iv in Figure 5C is oscillated vertically, fibers below the bead exhibit elevated I while fibers above the bead do not. In summary, the response of the fiber network to the oscillation of an embedded bead depends on the fiber network and the direction of oscillation.

DISCUSSION
We have shown a method to identify which fibers within a 3-D fiber network carry tension in response to the application of a localized force. Specifically, the method determines which fibers are oscillating in phase with an actively driven bead embedded within a pore of the network and assumes that pure elastic forces are driving these fiber motions. Here, fiber oscillation is measured by quantifying intensity fluctuations of all pixels in a plane confocal to the bead. As expected, Fourier magnitudes of pixel intensity fluctuations ( I ) at the drive frequency (F d ) are sensitive to even sub-pixel fiber oscillations ( Figure 1C). By looking at the I , we find that displacements are larger along the direction of the bead's oscillatory displacement (Figure 6), which is also the expected behavior in continuum viscoelastic materials [38]. In our system, however, we do not find quantitative agreement between the continuum-based theory and our data. Unlike in a continuum, we observe displacement propagation along complex paths at the scale of the mesh, likely dependent on the architecture of the local fiber network (Figures 5B,C). In addition to I , we find that pixels intensity fluctuations along fibers either have random phase (∠ I) or exhibit local pixel neighborhoods of similar ∠ I occurring on fibers that I indicate are under tension. In order to quantify and aid in the detection of such neighborhoods we compute the Orientational Order Parameter, or OOP. By use of the OOP, we identified distinctive sets of fibers that are engaged as the direction and amplitude of oscillation of applied force varied (Figures 4C,  5F,G). Importantly, while we are interested in the static, or zero-frequency force map, we opt to use low-frequency oscillation of the bead in order to take advantage of lock-in detection and further to provide resilience to forces beyond those applied by the optical tweezers (e.g., cell-generated forces). The choice of a low frequency ensures our findings are not complicated by nonlinear material properties of the system that emerge at higher frequency oscillations or viscoelastic properties of the local mesh. Our choice of applying small displacements relative to the scale of the pores also ensures  that we are measuring the current state of the system and not further complicate the measurement with additional network rearrangement and fiber alignment that emerge with larger applied deformations [21,22,39]. As compared to methods such as fiber tracing and application of static forces, our method may be more suitable for use in the presence of cells or other force centers, which can induce fiber displacements in addition to the actively driven displacements. Also, the lower frequencies and displacements are well within the technical capabilities of our optical tweezers-confocal microscope apparatus. Force transmission pathways as detected by fiber fluctuation analysis is highly variable from bead to bead due to the different network structure. One consistent observation is that fibers not appearing to be connected to the bead in the confocal 〈I〉 image still oscillated in phase with that bead (Figures 4, 5). This coordination is most likely due to the fact that not all force transmission pathways are in the plane containing the midpoint of the trapped bead. In order to identify potential out-of-plane force transmitting fibers, we overlay the OOP image on top of confocal z-stacks of the fiber network ( Figure 7). We find that there are fibers or a network of fibers which connect all the isolated in-plane high OOP and I regions to the bead. Black arrows in Figure 7 identify local sets of out-of-plane fibers connected on both sides to fibers exhibiting oscillations as detected by the local OOP (Figure 7). It is evident from the data that some fiber chains transmit forces while others do not. Extension of pixel fluctuation intensity analysis to include these out-of-plane fibers is a goal for future method development.
Lastly, we discuss a caveat of this intensity fluctuation analysis. I , ∠ I, and OOP should not be directly compared between fibers, especially fibers of different fluorescence brightness < I > . More precisely, a major factor determining the I is the gradient of fluorescence brightness along the direction of oscillation. For example, the gradients of 〈I〉 the along x ( Figure 8A) and y ( Figure 8B) axes for bead i in Figure 5 are shown. Magnified regions of Figures 8A,B show that a fiber aligned with the y axis has steep intensity gradients across the fiber in the x-direction and shallow gradient along the fiber in the y-direction ( Figures  8C,D). Consequently, detection of oscillation in the y-direction will not be as sensitive as in the x-direction. Of note, punctate bright regions ( Figures 8A,B black arrows) have steep intensity gradients along all directions and consequently our method should be equally sensitive to oscillations along any axis at such regions. We show that there are parts of the fiber network which have very shallow or non-existent gradient intensity. This lack of gradient is especially evident along a uniformly fluorescently labeled fiber, which makes it difficult to detect the displacement of that fiber. This observation motivates future work in which texture in fiber fluorescence intensity could be introduced to every fluorescent fiber thus mitigating anisotropic sensitivity to displacement. Such texture could be achieved either through a labeling strategy or the use of structured light illumination. A labeling strategy that introduces bright and high contrast images would enable the use of well-established tracking methods such as optical flow analysis, digital image correlation, or erosion-and-reconstruction methods [40][41][42][43][44][45][46][47][48][49].
In conclusion, we developed a straightforward approach to quantitatively determine which fibers are participants in force transmission highways. The use of local variations in the phase of intensity fluctuations provides sensitivity to even small displacements in the fiber network, where magnitude of intensity fluctuations were not distinguishable from noise. Ultimately, we aim to observe changes to the force transmission characteristics of biologically derived fiber networks and how these characteristics change as cells interact with the ECM. While we have shown the response of a fiber network in its relaxed state, future application of this technique lies in exploring how the response changes when the network contains cells. For example, cell contractile forces will apply a tensile preload that may change the subset of tensed fibers thus altering the scope of cell-cell force communication. Such observations of how force transmission pathways change in response to cellular contraction and ECM remodeling offers potentially new insights into processes such as growth, wound healing, cancer growth and metastasis.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation. Raw data related to the article can also be found at the following Dryad link: https://doi.org/10.7280/D1ZT21.

AUTHOR CONTRIBUTIONS
QH was responsible for planning, data acquisition, analysis, and writing for the paper. TM and AG were responsible for analysis and writing for the paper. AL and EB were responsible for planning, analysis, and writing for the paper.

ACKNOWLEDGMENTS
We would like to thank Micah Lim for assisting with data acquisition and editing the manuscript. We would like to thank Dr. Mark Keating and Alicja Jagiełło for their contribution in building the optical tweezers microscopy system. QH and EB acknowledges support from the United States Air Force Office of Scientific Research FA9550-17-1-0,193 and the Office of the President of the University of California. TM and AG acknowledges support from NIH R01 HL129008 and NIH T32 HL116270. AL acknowledges partial support from NSF-DMR-1709785.