Strain and rate-dependent neuronal injury in a 3D in vitro compression model of traumatic brain injury

In the United States over 1.7 million cases of traumatic brain injury are reported yearly, but predictive correlation of cellular injury to impact tissue strain is still lacking, particularly for neuronal injury resulting from compression. Given the prevalence of compressive deformations in most blunt head trauma, this information is critically important for the development of future mitigation and diagnosis strategies. Using a 3D in vitro neuronal compression model, we investigated the role of impact strain and strain rate on neuronal lifetime, viability, and pathomorphology. We find that strain magnitude and rate have profound, yet distinctively different effects on the injury pathology. While strain magnitude affects the time of neuronal death, strain rate influences the pathomorphology and extent of population injury. Cellular injury is not initiated through localized deformation of the cytoskeleton but rather driven by excess strain on the entire cell. Furthermore we find that, mechanoporation, one of the key pathological trigger mechanisms in stretch and shear neuronal injuries, was not observed under compression.

. Collagen hydrogel composition. v/v% of each component to produce a collagen gel concentration of 2.2 mg/mL. cylindrical gels of 1.25 mm height and 15 mm diameter. The center base of the mold was sterilized with 99% isopropyl alcohol, coated with a thin layer of vacuum grease, and covered with a 15 mm cover slip. Collagen-cell solution was inserted into the base cavity, sandwiched with a chemically activated, hydrophilic 25 mm top cover slip, and polymerized in a humidified cell incubator at 37 • C and 5% CO 2 for 45 minutes. Polymerized gels were removed with a custom-made tool designed to avoid gel damage during removal, and transferred to six-well plates filled with 3 mL of culture media. To maximize nutrient diffusion, the 15 mm cover slip was carefully removed before incubation. Collagen gel cultures were maintained in a humidified cell incubator at 37 • C and 5% CO 2 , and culture media was exchanged at 24 h and every 2 days thereafter. Experiments were initiated at day 7 in vitro (DIV).

Fluorescence Staining
For all compression experiments, neuron viability and geometry were assessed using calcein AM (Life Technologies), a membrane-permeant dye that is converted from non-fluorescent calcein AM to green-fluorescent calcein upon acetoxymethyl ester hydrolysis by intracellular esterases. Calcein AM has been shown to be well-retained within live cells, and has a high quantum yield. 2,3 Acute alterations in plasma membrane permeability following compressive loading were assessed by evaluating uptake of Alexa Fluor 568 hydrazide (AFH) (Life Technologies), a relatively low molecular weight (730 Da) and cell-impermeant fluorophore used for intracellular labeling experiments. 4 Used extracellularly, AFH cannot infiltrate the cell membrane until pores, approximately the size of AFH, develop due to mechanical insult. Previous studies have shown that this mechanoporation occurs even for significantly larger fluorescent dyes (MW 150 kDa) after mechanical stretch. 5 Prior to mechanical loading, cell cultures were placed into a Chamlide CMB chamber (Live Cell Instruments Inc., Seoul, South Korea), and incubated with calcein AM (2 µM) and AFH (10 µM diluted with dH 2 O) in prewarmed HA + 2% B-27 for 45 minutes. After rinsing with 1× PBS, a 15 mm cover slip was gently placed on top of the gel, and a solution of AFH (10 µM) in HA/B-27 was left in the chamber during mechanical insult and imaging. Studies have shown maintenance of neuronal cell viability and function at ambient CO 2 for up to 48 hours [6][7][8] in HA/B-27.

Confocal Microscopy and Live Cell Imaging
For cell compression experiments, three-dimensional image volumes of neurons before and immediately after mechanical insult were acquired using a Nikon A-1 confocal system mounted on a Ti-Eclipse inverted optical microscope controlled by Nikon NIS-Elements software (Nikon, Tokyo, Japan). To ensure physiological conditions, temperature within the microscope chamber was controlled at 37 • C using a closed-loop heater (Air-Therm heater, World Precision Instruments Inc., Sarasota, FL), and the high-speed compression device was allowed to thermally equilibrate for 2 -3 hours prior to imaging. A CFI Plan APO VC 60× water immersion objective (NA = 1.2; Nikon) was used for all mechanical compression experiments. Calcein AM-stained neurons were excited with an Argon ion laser (488 nm), and extracellular AFH was excited with a HeNe diode laser (561 nm). Confocal image stacks of 1024 × 1024 × 100 voxels (177 × 177 × 100 µm 3 ) were recorded with a z-step of 1 µm. Volumes were located at a depth of 100 -200 µm from the bottom of the collagen gel. For each experiment 6 -12 cells were imaged directly after mechanical insult, at four hours, and every hour thereafter until somatic lysis or leaking of calcein AM dye coupled with severe morphological changes occurred. Morphological changes and cell death were never observed before four hours in all experiments and thus imaging was begun then to avoid overexposure to the excitation laser.
Reflectance microscopy was used to visualize the fiber structure of the collagen gels. No stain is required and the samples were imaged with an APO 40× water immersion objective (NA = 1.15; Nikon) and Z-step of 0.63 µm, to yield an image volume of 1024 × 1024 × 176 voxels (212 × 212 × 110 µm 3 ).

Control experiments
Control experiments were conducted to confirm neuron population viability in the collagen hydrogels through 9 DIV (Fig.  S1). Culture preparation and maintenance remained the same as described earlier. Again, calcein AM was used as a live cell stain. Ethidium homodimer-1(EthD-1; Life Technologies), a cell-impermeant nucleic acid stain was used to stain dead cells. Upon compromise of the cell membrane, EthD-1 binds to nuclear DNA and, upon laser excitation, fluoresces red. Calcein AM-stained neurons were excited with an Argon ion laser (488 nm), and EthD-1 was excited with a HeNe diode laser (561 nm). Confocal image stacks of 512 × 512 × 200 voxels (635 × 635 × 248 µm 3 ) were recorded with a z-step of 1.24 µm. Neuron population viability was determined by manual counting aided by flNeuronTool. 9 Control experiments to assess the cell damage associated with the cell stains, and confocal microscopy imaging were performed under the same conditions as the impact experiments without applying compressive strain (Fig. S1).

3D Cell Compression Device
The 3D cell compression device applies mechanical deformation to neurons cultured in a 3D collagen gel. The primary goals of the cell compression device are (1) uniaxial compression across the 3D collagen cell culture with high strain rates and magnitudes that are characteristic of inertial TBI; [10][11][12][13] (2) displacement (as opposed to force) controlled deformation with well-characterized boundary conditions; (3) user defined temporal characteristics of the strain field; and (4) access to cells via confocal microscopy before, during, and after mechanical insult so that the immediate and long term biochemical events are measurable.
The cell compression device consists of a linear voice coil actuator (VCA) (LAS 13-18; BEI Kimco, San Marcos, CA) that is mounted on a custom frame attached to the confocal microscope stage (Fig. 1a). Uniaxial compressive strain is generated by linear motion of the VCA piston in contact with the top cover slip of the collagen gel (Fig. 1b). To ensure normal contact with the cover slip during deformation, the piston is fitted with a spherical tip. To overcome the limited range of motion of the piston and to increase sample accessibility, the VCA is mounted on a manual z-stage with a full range of motion of 100 mm. The sample chamber is situated in an independent x-y translation stage that is used to precisely center the sample with the piston. Prior to mechanical deformation, a custom cylindrical tool (with same diameter as the sample chamber) with a centered hemispherical cutout (with the same radius as piston tip) is inserted into the x-y stage. The VCA is lowered and the x-y stage is adjusted so that the piston tip fits into the centered cutout. Motion of the piston is produced via a closed loop control system consisting of a servo controller (SilverSterling S2-IG; QuickSilver Controls, Inc., Covina, California) and integrated hall-based displacement transducer. A computer is used to program and communicate with the controller for the desired piston motion. In the study, the VCA was programmed to deliver strain magnitudes of 0.30 and peak strain rates of 10 s −1 and 75 s −1 (Fig. 1e). Due to the inertia of the system, the actual measured strain differed slightly from the programmed.

Experimental Parameters
The uniaxial compressive strain field was validated under quasistatic loading on the 3D collagen neuron cultures (n = 7). Fluorescent microspheres (0.5 µm in diameter, carboxylate-modified, Life Technologies) were thoroughly mixed at 6 v/v% into the collagen gel during cell culture and maintained in the same conditions as the mechanical impact experiments. At day 7 in vitro, collagen cell cultures stained with calcein AM were mounted on the cell compression device and compressed by a 2% strain increment with a 2 minute loading time every 6 minutes to a total strain of 0.38. During the hold time, confocal image volumes of size 512 × 512 × 350 voxels (212 × 212 × 210 µm 3 ) were acquired within a depth of 100 -200 µm from the coverslip using an APO 40× water immersion objective (NA = 1.15, Nikon).

Full-field displacement computation
Full-field displacements of the collagen-embedded microspheres were tracked between loading steps using Fast Iterative Digital Volume Correlation (FIDVC). 14 FIDVC utilizes the iterative deformation method (IDM), which spatially refines the correlation window, allowing for accurate measurement of large, non-linear material displacements, which could arise near the neuron-collagen interface. The final interrogation window spacing, which sets the displacement spatial grid resolution was 8 voxels or 3.3 µm for measurements in x 1 − x 2 and 4.8 µm for measurements in x 3 . The displacement gradient tensor ∇ ∇ ∇u u u, computed with an optimal-9 tap finite difference kernel, was used to calculated the Lagrangian strain field E, given as The transverse to axial strain ratio ν was computed as the negative ratio of mean transverse strain E 11 and E 22 to axial strain E 33 averaged over all increments N, i.e.

Spherical Inclusion Problem
To determine the local strain values experienced by cells, a relation between far-field applied strain and cell-experienced strain must be understood. Cells compressed in a 3D matrix can be analytically compared (to first order) to an Eshelby inclusion solution given by Goodier, 15 which describes a spherical inclusion with radius a, Poisson's ratio ν I , and Young's modulus E I in an infinite, isotropic, linear-elastic solid matrix with Poisson's ratio ν M and Young's modulus E M subject to a uniaxial stress at infinity σ ∞ . The elastic mismatch between the matrix and inclusion gives rise to nonlinear displacements around the interface. The closed-form solution of the radial displacements u r and tangential displacements u θ solely due to the elastic mismatch is given by where A, B, and C are constants determined by equating the displacements and tractions at the interface, i.e. r = a. The values of the constants are where the matrix shear modulus µ M and inclusion shear modulus µ I are related to the Young's modulus and Poisson's ratio by The full-field displacements in the matrix are obtained by superposing the displacements in equation (3) with the solution for uniaxial compression/tension along x 3 , To compare between the experimentally determined displacement field around the cell and the Eshelby solution, an inclusion (a = 5 voxels, ν I = 0.49) in a matrix (ν M = 0.06) is subjected to a remote stress given by E ∞ 33 is the far-field strain applied by the cell compression device, and E ∞ 11 and E ∞ 22 are the strain components that arise due to the Poisson effect. The displacement along the loading axis u 3 is computed after transforming equations (3) -(6) to Cartesian coordinates for varying ratios of E I /E M . It is important to note that there are values of E I /E M that produce well-known deformations; E I /E M = 0 is the Eshelby solution of a spherical cavity subject to a remote stress; E I /E M = 1 is the solution for uniform strain; and E I /E M = ∞ is the Eshelby solution for a rigid inclusion subject to a remote stress. Qualitative comparison of the experimental local displacement field around neurons to the Eshelby inclusion solution for different ratios of E I /E M showed good agreement within a range of 0.5 to 2.0 (Fig. S2), suggesting that the imposed far-field strain is not significantly affected by a local elastic mismatch, which is consistent with research on other cell types. 16 Thus, cells are assumed to locally undergo the imposed far-field strain during impact.

Image Preprocessing
Converting raw confocal images of neurons to quantitative metrics that can be used for injury analysis requires a multiphase image processing routine (Fig. S3). The 3D confocal micrographs were split to isolate the calcein AM and Alexa Fluor Hydrazide channels independently and were converted to .tif files, which were then loaded into MATLAB (Natick, MA) for preprocessing. Images were rigidly translated to correct for experimental drift by maximizing the cross-correlation between an image and the image at the next time point using the calcein AM channel.
A median filter with a 3 × 3 × 3 voxel window size was then applied twice for noise reduction, followed by a gamma correction with a coefficient of 0.7 (I(x x x) = I(x x x) 0.7 ) for contrast enhancement. To remove background noise, all intensity values below the mode of the image were set to zero.

4/14
Building Neuron Tree Structure The tree data structure, commonly used in graph theory, is a powerful way to describe the morphology of a neuron. In this study, the neuron tree structure, T = {V, E}, is composed of a set of n nodes or vertices V connected by a set of n − 1 edges E that are directed away from a root that has no parent node (Fig. S4a). Nodes are classified by the degree, d, or number of children that directly follow, i.e. d = 1, extension; d > 1, fork; and d = 0, end or leaf. For quantitative geometric and injury analysis, nodes and edges can store important information (e.g. position, radius, and strain). Implementation of the neuron tree allows for easy and efficient traversal along paths by eliminating repetition.
To build a tree structure, each calcein AM image was exported as .tiff format and loaded into Neuron Studio, a software that performs semi-automatic tracing and reconstruction of neuron structures from confocal images. 17 From a manually selected seed point at the soma center, Neuron Studio automatically thresholds and skeletonizes the image with a topology and geometry-preserving algorithm that computes the medial axis. 18,19 Branch diameters are estimated at each node along the skeleton or medial axis using the Rayburst sampling algorithm. 18,19 In this algorithm, the image is trilinearly interpolated along ray vectors that are cast from a node, and are allowed to grow until a specified intensity is reached. This intensity value is automatically adjusted by Neuron Studio according to the local intensity distribution of the image. Final node diameters are determined via rayburst spheres, estimated by computing the median ray vector distances for each node, which are less affected by local surface irregularities.
The Neuron Studio tracing and reconstruction procedure produced a tree structure of the neuron (Fig. S4b), which was carefully edited for validity, i.e. extensions automatically built from other neurons in the image were carefully removed using the built-in 3D viewer. To prevent undersampling, the "Discretization Ratio" parameter, which sets the sampling density of nodes, was set to the minimum value of zero, producing the most dense nodal connectivity. Additionally, spurious extensions that were produced from image noise were removed by setting the "Minimum Length" parameter to 5 µm. The image channel of the cell-impermeant Alexa Fluor Hydrazide dye was manually traced since the neuronal extensions were not as well defined. In most cases, the tree consisted of a soma and one to two short neurites.
After editing, the neuron trees were exported as 'swc' format, read into MATLAB, and converted to a tree structure for quantitative injury analysis. Functions that were not custom written for a specific application came from the TREES toolbox, 20 a graph theory library specifically for neuron tree structures in MATLAB.

Defining local coordinate frames for neuron structures
Defining quantitative injury metrics requires the definition of a local coordinate frame. Each neurite extension, if treated like a smooth parameterized space curve, x x x(s), can be described by an orthonormal Frenet frame, f f f (s). 21 At each point on the curve, the Frenet frame is defined as where t t t(s), n n n(s), and b b b(s) are the tangent, normal, and binormal unit vectors (Fig. S5a). They are mathematically expressed as The curvature, κ(s), and twist, τ(s), of x x x(s) are related to t t t(s), n n n(s), and b b b(s), and their derivatives by the Frenet-Serret system of differential equations, but are not utilized in this study as they are dependent on the curve fit third derivative x x x (s), which is not guaranteed to be smooth for our chosen cubic spline fitting. Although the Frenet frame f f f (s) can be computed along the extensions of the cell, its rotation about t t t(s) can behave erratically at degenerate points where the curvature vanishes, e.g. inflection points or straight segments. To circumvent this issue for injury analysis, the injury metrics proposed in this study will be invariant to rotation about t t t(s).

Practical Implementation of a discrete local coordinate frame for neuron structures
In practice, using equations (8) to define a local coordinate frame is numerically challenging for many reasons. Extensions are not made of smooth curves, but piecewise linear functions that can sharply change direction at certain nodes due to image noise. Additionally, node sampling is highly dependent on the tracing parameters in Neuron Studio. To circumvent these issues, a series of processing steps was performed before computing the local coordinate frame.

5/14
Tree resampling: To begin the processing procedure, the neuron tree was resampled using the 'resample tree' function in the TREES toolbox, which redistributes the nodes at a user defined equal distance away from each other. The resampling distance was chosen as the 25th percentile of all edge lengths, such that most geometric features are preserved. Tree smoothing: After resampling, trees were smoothed to eliminate sharp corners using Taubin smoothing, an algorithm commonly used to remove noise on polygonal meshes while preserving salient shape features. [22][23][24] To adapt this algorithm for trees, first the discrete Laplacian of x x x is computed for each extension node (d = 1) with the expression Geometrically, ∆x x x i are the average positions of neighboring nodes to the ith node. In general there are more advanced weighting schemes for the discrete Laplacian operator. However, since the tree is resampled at equidistant points, the effects would be minimal. A Laplacian filter is then applied by updating the position x x x with ∆x x x, given by the equation where λ is a scale factor that controls the extent of smoothing, and is 0 < λ < 1. The Laplacian filter must be applied iteratively many times to produce significant smoothing. However, doing so would result in large tree contraction or shrinking. To overcome this issue, Taubin incorporated an expansion step following each contraction step. Mathematically, the expansion step is similar to equation (11), but λ is replaced with a negative scale factor µ such that µ < −λ . The resulting transfer function of the expansion and shrinking operation has a pass-band frequency given by All trees in this study were smoothed 100 times with λ = 0.63 and k PB = 0.1, chosen from. 24 Computing discrete tangent: Having sufficiently smooth trees allows for computation of f f f (s), for which there are many numerical methods. 25,26 In this study, a piecewise cubic spline (C 2 differentiable) was fit to a set of tree paths using the "spline" function in MATLAB. Each path was defined as set of extension type nodes that connect a root and fork, root and end, fork and fork, or fork and end. As opposed to other methods, cubic splines return a set of third order polynomial coefficients that can be analytically evaluated to compute the geometric features at any point.

Image Post-Processing
To accurately tag the formation of blebs or to determine when cell death occurred, all image information outside of the tree was removed (Fig. S6a). To perform this masking procedure, connected nodes in the tree were isolated and converted to a nodal mask by setting interior voxels to 0 and exterior voxels to 1 (Fig. S6 b and c). For each nodal mask, a convex hull (CH) was computed and superimposed with the previous nodal mask using an OR logical operation to construct the CH of the whole cell (Fig. S6d). The final image mask was computed by dilating the CH of the whole cell using a Euclidean distance transform, the 'bwdist' function, where all values within a distance of 32 voxels were set to 0 and outside were set to 1 (Fig. S6e).

Local Strain Transformation
The strain that neurons experience is highly dependent on the local orientations of their neurites. For example, a neurite that is aligned with the principal loading axis will experience more axial strain than one that is transversely aligned with the loading axis. Analysis for this orientation-dependent deformation draws from classic laminate theory and strain rosette calculations, where the global strain or stress is transformed to coordinates of the material fibers or strain gauges. For both of these cases, the analysis is simplified using plane stress assumptions, which is not valid for neurons. To extend the theory into 3D, the peak far-field strain E ∞ must be rotated to the direction of each edge in the neuron tree. The transformation is expressed as where E ( f f f ) is strain matrix in the local basis, Q is the rotation matrix. E ∞ is given from the strain applied by the cell compression device and the transverse to axial strain ratio ν as Using the Frenet frame ( f f f ) along x x x(s), Q is given by · e e e 1 b b b · e e e 2 b b b · e e e 3 n n n · e e e 1 n n n · e e e 2 n n n · e e e 3 t t t · e e e 1 t t t · e e e 2 t t t · e e e 3  

6/14
where e e e is the global Cartesian coordinate frame, t t t is the local tangent vector, and e e e 3 is the loading axis. Applying the rotation in equation (15) to E ( f f f ) , equation (13) becomes t t t n n n · E ∞ · n n n n n n · E ∞ · t t t t t t · E ∞ · t t t.
Since t t t is the most uniquely defined geometric quantity in the system, we compute the axial strain, defined as ( Fig. S5 d and e). Because the orthonormal vectors n n n and b b b are non-unique in our system, we define an average axial shear strain as After equations (17) and (18) are computed for every point along a neuron, the mean axial strain is defined as and the mean axial shear strain is defined as where L is the total length of the neuron, i.e. L = x x x(s) ds.

Injury Assessment
Bleb Formation: For each neuron that exhibited clear blebbing, a time series of maximum intensity projection images was created. The image at the time point at which severe blebbing occurred was loaded into flNeuronTool, a software for 3D visualization and proofreading of neuronal morphologies. 9 All blebs were tagged using the 'probing' mode and imported into MATLAB for analysis. The maximum compressive axial strain and axial shear strain of all nodes within a specified radial distance of each bleb were computed. This distance was equal to 3.35 ± 0.73µm, computed as the mean bleb radius of 26 blebs as manually recorded from three representative neurons.
Membrane Poration: AFH trees were loaded into MATLAB and each node sphere along the neuron tree was converted into a binary mask. Neurons were considered porated if the mean intensity within the mask increased to the mean background intensity. The same analysis was performed for the membrane poration controls.
Cell Death Assessment: Time of cell death was assessed as the timepoint at which somatic lysis or a combination of leaking of live cell dye calcein AM and severe morphological deformities occurred. Leaking of calcein AM fluorescence has been shown to correlate with other cell death markers. 27 Bivariate Compressive Axial Strain Fitting and Probability Density Estimation: The scatter data of the mean compressive axial strain E c as a function of time at cell death t d was fit to a bivariate normal distribution using the minimum covariance determinant (MCD) estimator 28,29 implemented in the LIBRA MATLAB library. 30 MCD is commonly used in bivariate fitting since outliers in the data can greatly influence the fitting. The major and minor axis of the 95 % ellipse were computed from the eigenvalues of the covariance matrix of the fit. To qualitatively confirm the validity of the bivariate normal distribution, the 2D probability density was computed using the commonly used non-parametric adaptive algorithm for kernel density estimation. 31