Quantification of Cardiomyocyte Alignment from Three-Dimensional (3D) Confocal Microscopy of Engineered Tissue

Abstract

Biological tissues have complex, three-dimensional (3D) organizations of cells and matrix factors that provide the architecture necessary to meet morphogenic and functional demands. Disordered cell alignment is associated with congenital heart disease, cardiomyopathy, and neurodegenerative diseases and repairing or replacing these tissues using engineered constructs may improve regenerative capacity. However, optimizing cell alignment within engineered tissues requires quantitative 3D data on cell orientations and both efficient and validated processing algorithms. We developed an automated method to measure local 3D orientations based on structure tensor analysis and incorporated an adaptive subregion size to account for multiple scales. Our method calculates the statistical concentration parameter, κ, to quantify alignment, as well as the traditional orientational order parameter. We validated our method using synthetic images and accurately measured principal axis and concentration. We then applied our method to confocal stacks of cleared, whole-mount engineered cardiac tissues generated from human-induced pluripotent stem cells or embryonic chick cardiac cells and quantified cardiomyocyte alignment. We found significant differences in alignment based on cellular composition and tissue geometry. These results from our synthetic images and confocal data demonstrate the efficiency and accuracy of our method to measure alignment in 3D tissues.

INTRODUCTION

Microstructural organization is essential to the organogenesis and function of tissues and organs. The alignment of cells and matrix components along a preferential axis contributes to the acquisition of mature geometries, mechanical integrity, electrical conduction, and contractile efficiency. The importance of alignment has been studied in various tissues including bone (Banglmaier et al., 2015), ocular tissue (Avila & Bueno, 2015), neurite growth (Singh et al., 2016), vasculature (Gasser et al., 2006), and the heart (Streeter et al., 1969; Gilbert et al., 2007). Tissues adapt to environmental stimuli by remodeling and reorganizing their microstructure to optimize performance. Cardiomyocytes (CM) cultured in vitro align parallel to the stretch direction (Matsuda et al., 2005; Nguyen et al., 2013) and endothelial cells (EC) elongate and align in the direction of steady flow or cyclic stretch (Galbraith et al., 1998; Kaunas et al., 2005). Genetic errors associated with abnormal myofiber alignment and fibrosis form the basis of a range of clinically important cardiomyopathies (Carrier et al., 1997; Akdis et al., 2016; Finsterer & Stollberger, 2016) and are associated with abnormal myocardial structure and function in congenital heart disease (Hahn et al., 2015).

The architecture of the heart has been widely studied, detailing a laminar organization of helically arranged CM with similarly aligned support cells, collagen, and other extracellular matrix components (Gilbert et al., 2007). Establishing this myoarchitecture during normal cardiac morphogenesis requires complex temporal and spatial myocardial maturation (Sedmera et al., 2000; Tobita et al., 2005). The alignment and organization of cardiac tissue is a primary factor in generating the stress and strain required during contraction and relaxation (Waldman et al., 1988; Omens et al., 1991; Costa et al., 1999; Ashikaga et al., 2009). Propagation of action potentials also requires alignment of adjacent CM, as conduction velocities are several times faster in the longitudinal versus transverse direction (Panfilov & Keener, 1993; Punske et al., 2005; Hooks et al., 2007; Dhein et al., 2014). Structural organization is established during early embryonic development, where myofiber orientations proceed through several iterations (Tobita et al., 2005). When the embryonic heart is exposed to abnormal mechanical loading, myofiber organization is disrupted, indicating biomechanical regulation (Sedmera et al., 2000; Tobita et al., 2005) and producing altered developmental trajectories (Burggren & Reyna, 2011). In the adult, various disease states can disrupt cardiac structural organization. The border zone in myocardial infarction (MI) exhibits myocardial disarray, which can lead to reentrant arrhythmias, slow conduction, and electrical instability (Wickline et al., 1992; Sosnovik et al., 2009; Rutherford et al., 2012). In the hypertrophic heart, the transmural helical variation of CM is reduced and tissue organization is diminished (Roberts & Ferrans, 1975; Schmitt et al., 2009; Hill et al., 2014). Myocardial fibrosis in patients with hypertrophic cardiomyopathy is associated with increased electrical heterogeneity and, importantly, significant risk of sudden, arrhythmic death (Kadish et al., 1988; Valderrabano et al., 2001; Ripplinger et al., 2007; Muthappan & Calkins, 2008).

CM loss after MI impairs cardiac function, leading to a 50% 5-year survival rate for heart failure patients, in part because the heart lacks the capacity to restore lost CM (Mozaffarian et al., 2016). Implantable cells and engineered cardiac tissues (ECTs) for cardiac repair have shown positive preclinical and clinical results and represent promising therapies for heart failure (Bolli et al., 2011; Sanganalmath & Bolli, 2013; Masumoto et al., 2014, 2016; Fisher et al., 2015). Cardiac tissue engineering approaches attempt to recapitulate the process of CM organization to generate constructs with structural and functional properties optimal for implantation. Substrate stiffness (Ribeiro et al., 2015; Li et al., 2016), micropatterned coatings (Feinberg et al., 2012; Tsang et al., 2015), and tissue geometry (Bian et al., 2014; Nguyen et al., 2014) all influence CM size, shape, and alignment. In addition, mechanical and electrical conditioning may accelerate tissue organization (Nunes et al., 2013; Stoppel et al., 2016). Quantifying alignment in engineered tissues is therefore critical to understand the maturation and viability of these constructs and optimize culture methods.

There are a variety of methods for quantifying alignment in two-dimensional (2D) images. Measures of preferential cellular orientation are often determined from 2D histologic images. One of the earliest and still widely used is the Fourier transform radial sum, which calculates a single, global 2D alignment measurement (Chaudhuri et al., 1987; Petroll et al., 1993). Fitting an ellipse model to cell boundaries can quantify individual orientations, but is highly dependent on the segmentation method (Bray et al., 2010; Xu et al., 2011; Nectow et al., 2014). In cases where ellipse models are not appropriate, segmented images have been processed using the Hough transform (Bayan et al., 2009) or eigenanalysis of the 2nd-order moment tensor (Vader et al., 2009) to determine local orientations. In these methods, the orientation operation is performed on image subwindows, where the size of the subwindow defines “local.” Methods that do not require segmentation rely on the image intensity gradient, where the gradient magnitude is large at object boundaries and its direction is normal to the boundary. Subwindow analysis of image gradients has been demonstrated as an accurate method for determining local orientations (Chaudhuri et al., 1993; Hong et al., 1998; Karlon et al., 1999; Grosberg et al., 2011; Sun et al., 2015). Our group has applied this method to quantify myofiber alignment in the embryonic heart and in ECTs (Tobita et al., 2005; Masumoto et al., 2016). In all of these methods, the distribution of local orientations is used to determine alignment, which measures the amount that structures within an image are oriented parallel to a given direction (usually the mean direction).

There are fewer methods to quantify cell or matrix orientations in three-dimensional (3D) tissues, given the increased complexity, time required, and cost for 3D image data acquisition and analysis. Imaging modalities such as diffusion tensor magnetic resonance imaging (DTMRI; Hsu et al., 1998) and polarized light microscopy (Tower & Tranquillo, 2001) provide direct information on tissue orientation but cannot distinguish cell types. The resolution of DTMRI may not be sufficient for small scale tissues. Segmentation of cells or fibers can provide 3D orientation data, but 3D segmentation is costly and subjective. An automated tracking algorithm developed for analysis of segmented collagen images (FIRE; Wu et al., 2003) has been combined with the curvelet transform to first remove noise from the image and then quantify tissue features including fiber width, length, density, and orientation (CurveAlign; Bredfeldt et al., 2014 a , 2014 b). Segmentation quality and the >20 required parameters limit the utility of this method. Eignenanalysis of the Hessian matrix calculates orientations at the voxel level and does not require a segmented image (Frangi et al., 1998; Daniels et al., 2006). This method is improved when applying an adaptive filter size, though it does significantly increase computation time (Daniels et al., 2006). Other methods rotate local subregions to maximize the z-axis 2nd-order moment (Reuze et al., 1993), perform convolutions with anisotropic Gaussian filters aligned in different directions, where the orientation is given by the direction that produces the maximal result (Robb et al., 2007; Wirjadi et al., 2009), or compute the distances to the fiber boundary in all 26 directions to identify the major inertia axis (Altendorf & Jeulin, 2009; Altendorf et al., 2012). These methods have varying degrees of accuracy and require some additional validation.

Methods for 3D orientation based on the image intensity gradient are advantageous as they avoid segmentation and the gradient is faster to compute than Hessian matrices. However, these approaches are sensitive to noise, the method used to compute the gradient, and subregion size (if one is used), and therefore require validation using synthetic data. The structure tensor (T) is the tensor product of the gradient vectors (Weickert, 1998; Puspoki et al., 2016). T is usually summed or averaged (by convolution with a Gaussian) over a region to provide a localized covariance matrix of the gradient. Eigenanalysis of T provides information about the shape and orientation of the local region. The structure tensor has been used to quantify myofiber orientations in serial reconstructed sheep atria images (Zhao et al., 2012 a, 2012 b) and in high-resolution 3D cardiovascular magnetic resonance images of explanted rat hearts (Gilbert et al., 2012; Bernus et al., 2015). These methods identified laminar orientation well, but may not be applicable to myocyte orientation due to resolution limits. In orthopedic studies, T has been used to quantify anisotropy in trabecular bone (Tabor & Rokita, 2007; Kersh et al., 2013; Larsson et al., 2014). These studies used quantitative computed tomography (CT) and magnetic resonance imaging and compared well with the gold-standard mean intercept length method. However, the results were less accurate for clinical CT resolution (Tabor & Rokita, 2007; Kersh et al., 2013; Larsson et al., 2014).

In order to investigate cell and myofiber alignment within 3D ECTs, we applied structure tensor analysis to 3D confocal images and demonstrate an accurate method to compute local orientation of tissue structures and quantify alignment. Our method incorporates an adaptive subregion technique to account for different size scales. For validation, we tested our orientation measurements and alignment metrics on synthetic images. We then analyzed whole-mount ECTs to measure CM orientations identified using cardiac troponin T immunostaining. We applied tissue clearing to enable nondestructive imaging of thick tissues and defined a local coordinate system based on whole-tissue geometry. Our method provided robust measurement of 3D CM alignment to enable the assessment of the impact of ECT geometry, cell source, and electrical conditioning on tissue alignment. Our results suggest this method is applicable to other tissues.

MATERIALS AND METHODS

Coordinate System

The spherical coordinate system used in this study defines azimuthal angle θ as the counterclockwise angle from the positive x-axis to the orthogonal projection on the xy-plane and polar angle φ measured from the positive z-axis. Spherical coordinates are given as (θ,φ).

Local 3D Orientation

Our method for determining local orientations of objects in a 3D histology image, I, relies on the image intensity gradient, G=[G _x ,G _y ,G _z ] ^T , at each voxel (equation 1). We computed G based on convolution with Gaussian kernels (h, equation 2), similar to Karlon et al. (1999), where σ controls the area of influence (5 μm for all results in this study) and the kernel size (s, −s≤x,y,z≤s) was chosen such that the minimum value of h was 0.01.

(1) $$G_{x} {\equals}h_{x}\, {\ast}\,I,\,G_{y} {\equals}h_{y} \,{\ast}\,I,\,G_{z} {\equals}h_{z} \,{\ast}\,I,$$

(2) $$\!h_{i} \left( {x,y,z} \right){\equals}{{2i} \over {\sigma ^{2} }}{\rm exp}\left( {{{{\minus}\left( {x^{2} {\plus}y^{2} {\plus}z^{2} } \right)} \mathord{\left/ {\vphantom {{{\minus}\left( {x^{2} {\plus}y^{2} {\plus}z^{2} } \right)} {\sigma ^{2} }}} \right. \kern-\nulldelimiterspace} {\sigma ^{2} }}} \right),\,i\in\left\{ {x,y,z} \right\}.$$

We then computed the local dominant orientation within image subregions of size m×m×m. For each subregion, we collected all of the gradient vectors as the set G ^m =[G ₁, … , G _n ] and computed the structure tensor T (equation 3). As we weighted the gradient vectors equally, we simply summed over the subregion. We computed the eigenvalues (|λ ₁|≤|λ ₂|≤|λ ₃|) and corresponding eigenvectors (e ₁, e ₂, e ₃) of T. The eigenvector e ₁, with the smallest magnitude eigenvalue, is the local orientation. This analysis treats G ^m as a girdle distribution in spherical statistics. A girdle distribution of axes is concentrated around a great circle; it is an equatorial distribution (Fisher et al., 1987). The polar axis of a girdle distribution is the vector normal to the equatorial plane. As image gradient vectors are normal to the surfaces of objects within the image, this polar axis is thus the dominant orientation of objects within the local subregion. The polar axis is the eigenvector associated with the smallest eigenvalue of the structure tensor T (Fisher et al., 1987).

(3) $${\bf{T{\equals}G}} ^{m} {\bf{G}} ^{{m^{T} }} .$$

Subregion size m defines “local” and is an important parameter in our method. Since biological objects (cells, nuclei, cytoskeletal elements, collagen fibers, etc.) can have a range of sizes, even within a single image, using a fixed m is not ideal. Therefore, we applied an adaptive subregion scale, based on the eigenanalysis of T. We first set minimum and maximum subregion sizes (m _min and m _max, respectively) and divided the image into evenly spaced subregions of size m _min. For each subregion, keeping the centroid fixed, we computed the local T from G ^m , as described above, for m values between m _min and m _max. At each size m, we calculated the shape parameter (γ, equation 4) based on the eigenvalues of T. We selected the subregion size m that produced the minimum γ. Small γ values indicate blob or rod-shaped objects (likely CM) as opposed to plate-like objects (Fisher et al., 1987). For all of the results in this study, we used an m _min of 30 μm and m _max equal to one-fourth the smallest image dimension. We increased m in 10 μm steps.

(4) $$\gamma {\equals}{{{\rm In}\left( {{{\lambda _{3} } \mathord{\left/ {\vphantom {{\lambda _{3} } {\lambda _{2} }}} \right. \kern-\nulldelimiterspace} {\lambda _{2} }}} \right)} \over {{\rm In}\left( {{{\lambda _{2} } \mathord{\left/ {\vphantom {{\lambda _{2} } {\lambda _{1} }}} \right. \kern-\nulldelimiterspace} {\lambda _{1} }}} \right)}}.$$

As an additional step, after dividing the image into subregions of size m _min, we rejected any subregion with an insufficient mean intensity. This step removes subregions that contain no objects. For this study, we set the intensity threshold based on the global mean intensity of the entire 3D image. For our validation tests, we set the threshold at 1.5 times the global mean intensity and for our ECT images we set the threshold at five times the global mean intensity. Our method requires five input parameters: (1) the standard deviation of the gradient kernel, (2) m _min, (3) m _max, (4) the subregion step size, and (5) the intensity threshold. The major steps of our method are outlined in Figure 1.

Figure 1 Flow chart depicting the major alignment method steps. After preprocessing, the image intensity gradient is computed. The image is divided into subregions of the set minimum window size (m _min). At each subregion centroid, eigenanalysis of the structure tensor for window size m gives the polar axis and shape parameter γ. This process is repeated for m between m _min and m _max. The polar axis with the minimum γ is defined as the local orientation. After computing the local orientation for each subregion, κ and OOP are calculated for the image.

3D Alignment in Whole-Mount Histology Images

The local orientations computed above form the set V=[v ₁, … , v _n ], where n is the total number of orientations. Several methods for measuring alignment in histology images have been proposed. These metrics attempt to determine the degree to which objects are oriented along a common direction. In 2D circular statistics, angular standard deviation can be used to define alignment, where larger values indicate a more isotropic organization (Fisher, 1995). The most common measurement in 2D, however, is the orientational order parameter (OOP, equation 5; Drew et al., 2015). An OOP of −1 indicates perfect alignment perpendicular to the mean direction, 0 is random alignment, and 1 is perfect alignment parallel to the mean direction. OOP is based on the cosines between the mean direction and the local orientations (ψ).

(5) $${\rm OOP}{\equals}{2 \over n}\mathop{\sum}{{\rm cos}^{2} } \left( \iPsi \right){\minus}1.$$

Spherical statistics provides an alignment measurement for 3D samples. The Watson bipolar distribution describes axial data, where the positive and negative directions of a vector are equivalent. Bipolar data is predominantly oriented along a principal axis. The concentration parameter (κ) measures the dispersion around the principal axis; larger κ indicates greater alignment (Fisher et al., 1987). We treated our set of dominant local orientations, V, as a Watson bipolar distribution and computed κ as our alignment measurement (equation 6). The calculation of κ is based on the largest eigenvalue (λ ₃) of the structure tensor computed from V. We additionally computed the principal axis (μ), which is the eigenvector associated with the largest eigenvalue (Fisher et al., 1987).

(6) $${{\lambda _{3} } \mathord{\left/ {\vphantom {{\lambda _{3} } n}} \right. \kern-\nulldelimiterspace} n}{\equals}{{\mathop{\int}\limits_0^1 {x^{2} {\rm exp}\left( {kx^{2} } \right)dx} } \!\mathord{\left/\! {\vphantom {{\mathop{\int}\limits_0^1 {x^{2} {\rm exp}\left( {kx^{2} } \right)dx} } {\mathop{\int}\limits_0^1 {x^{2} {\rm exp}\left( {kx^{2} } \right)dx} }}} \right. \kern-\nulldelimiterspace} {\mathop{\int}\limits_0^1 {x^{2} {\rm exp}\left( {kx^{2} } \right)dx} }}.$$

We compared κ with the OOP computed for a set of test images. To compute the OOP, we used cos(ψ _i )=v _i ·μ. In addition, we computed the 95% confidence cone major angle (α) surrounding the principal axis. To compute the confidence cone, we followed the procedure detailed by Fisher et al. (1987, §6.3.1(ii). As alignment increases, α becomes smaller. The 95% confidence cone is sensitive to the number of orientation vectors. To keep the number of vectors equal between samples, we randomly sampled (with repetition) 100 vectors from V and computed α. We repeated the sampling 25 times and calculated the mean α.

Validation

We tested our 3D orientation method using synthetic fiber images that included added noise to validate μ and κ measurements. Each image was comprised of 200 cylindrical fibers embedded in a 512×512×200 voxel domain, where the voxel size was (1.24×1.24×3.0 μm). Fiber lengths and diameters were selected from normal distributions with a mean of 100 and SD of 20 μm for lengths and a mean of 15 and SD of 2 μm for diameters. We selected the mean fiber diameter based on measurements from our confocal data. Fiber centerline axes were sampled from a Watson bipolar distribution (Chen et al., 2015). We varied the principal axis (μ ₀) or concentration parameter (κ ₀) of the sampling Watson distribution to test the accuracy of our method. We expect that the local orientations measured by our method have a μ and κ similar to μ ₀ and κ ₀. To test μ, κ ₀ was fixed to 3.0 and we varied the azimuthal (θ) and polar (φ) angles of μ ₀ between −180° and 180° and 10° and 90°, respectively, and generated one image for each μ ₀ (25 total images). To test κ, μ ₀ was fixed to (135°, 45°) and we varied κ ₀ between 0.25 and 10 (27 total images). We added random, normally distributed noise (mean 0, SD 2 μm) to the centerlines in order to incorporate fiber “waviness.” Fibers were added as solid objects to create a binary image, to which we then added white Gaussian noise (mean 0, variance 0.02). Before computing local orientations, we applied 3D Gaussian blur (SD of 3 μm) to replicate the preprocessing step from the confocal image analyses.

We computed local orientations for each test image using our algorithm described above. To validate μ by our method, we tested for a common principal axis between the orientations measured by our method (V) and the set of prescribed fiber centerline axes (V ₀). We performed this test for each image (n=25) and considered a probability value (p)<0.05 significant to reject the hypothesis of a common principal axis (Fisher et al., 1987). We then examined correlation between the set of measured principal axes (M=[μ ₁, … , μ ₂₅]) and prescribed principal axes (M ₀=[μ _0,1, … , μ _0,25]), where a coefficient (ρ) of 1 indicates perfect correlation (Fisher et al., 1987). We also computed the R ² value for the fit of the identity line (y=x) to the measured versus prescribed μ components. To validate κ, we tested whether the sets of V and V ₀ had equal concentrations (Fisher, 1986). This test was performed for each image (n=27) and we considered p<0.05 significant to reject equal κ. We then calculated the Pearson’s coefficient (r) to examine linear correlation between the 27 measured κ and prescribed κ ₀ values. In addition, we computed the R ² value for the fit of the identity line to the measured versus prescribed κ values.

One goal of our method is to quantify the degree of myofiber alignment within 3D histologic images from various formulations or culture conditions. Therefore, we tested the robustness of our alignment metric, κ. We selected six values of κ ₀ for the sampling Watson distribution (1, 2, 3, 6, 8, and 10) and generated five noisy fiber images for each one. For this test, we only included 50 fibers per image to reduce the time cost. The standard deviations of the normal distributions used to select fiber lengths and diameters were both 0. To demonstrate that our alignment measurements were independent of orientation, we varied μ ₀ of the sampling Watson distributions so that the five images had different principal axes. We computed local orientations and κ for each image using our method. For each group of five images with the same input κ ₀, we tested if κ for all of the measured orientation distributions were the same (Fisher, 1986). A p<0.05 rejected the hypothesis of a common κ among all five images.

ECT Generation

Our methods for fabricating linear ECTs using embryonic chick (Tobita et al., 2006), embryonic rat (Fujimoto et al., 2011), and human-induced pluripotent stem cells (h-iPSC)-derived cardiac cells (Masumoto et al., 2016) have been published. In brief, we generated ECTs from h-iPSCs that have been induced to differentiate into CM, EC, and mural cells (Masumoto et al., 2016). Embryonic chick CM were harvested from day 10 white Leghorn chick embryos. Collected cells (~3.0×10⁶ cells/ECT) were mixed with acid-soluble rat-tail collagen type I (Sigma, St. Louis, MO, USA) and matrix factors (Matrigel; BD Biosciences, Franklin Lakes, NJ, USA). Cell/matrix mixture was performed as follows. (1) Cells were suspended within a culture medium [high (25 mM) glucose modified Dulbecco’s essential medium; Life Technologies, Carlsbad, CA, USA) containing 20% fetal bovine serum (Life Technologies). (2) Acid-soluble collagen type I solution (pH 3) was neutralized with alkali buffer (0.2 M NaHCO₃, 0.2 M HEPES, and 0.1 M NaOH) on ice. (3) Matrigel (15% of total volume) was added to the neutralized collagen solution. (4) Cell suspension and matrix solution were mixed with a final collagen type I concentration of 0.67 mg/mL. Linear embryonic chick or h-iPSC ECTs were constructed using a collagen type I-coated silicone membrane six-well culture plate (TissueTrain; Flexcell International, Hillsborough, FL, USA) and FX-5000TT system (Flexcell International). The center of the silicone membrane of a TissueTrain culture plate was deformed by vacuum pressure to form a 20×2 mm trough using a cylindrical loading post (FX-5000TT); ~200 μL of cell/matrix mixture was then poured into the trough and allowed to gel for 120 min in a standard CO₂ incubator (37°C, 5% CO₂) to form a cylindrical construct. Nylon mesh anchors attached to the TissueTrain culture plate held both ends of the construct in place. Once the tissue gelled, the culture plate was filled with preculture medium [α minimum essential medium (Life Technologies) supplemented with 10% fetal bovine serum (FBS), 5×10⁻⁵ M 2-mercaptoethanol (Sigma) and 100 U/mL Penicillin-Streptomycin (Life Technologies)]. Constructed ECTs were cultured for 14 days with medium change every other day (Fig. 2a).

Figure 2 Example engineered cardiac tissue (ECT) geometries and whole-mount images. a: Linear ECT geometry constrained at the long axis ends by mesh nylon anchors, forming a cylindrical shape. b: 40× paraffin section of a day 14 linear ECT. Cardiomyocytes (CM) are predominantly located at the outer surface. c: Projection of a 10× three-dimensional (3D) linear ECT whole-mount confocal image. d: Large-format ECT (LF-ECT) formed by a polydimethylsiloxane (PDMS) mold. e: LF-ECT removed from the mold. Thin, cylindrical bundles connect at junctions. Green box shows a typical imaging location. f: 40× image of a whole-mount LF-ECT acquired at the bundle midsection. The specificity of 4',6-diamidino-2-phenylindole (DAPI) and cardiac troponin T (cTnT) staining is similar to paraffin sections (b). g: Projection of a 10× 3D LF-ECT whole-mount image. h–j: xy (h), yz (i), and xz (j) views of the LF-ECT in (g). Dashed lines show locations of each section. k, l: Attenuation correction of the cTnT image channel showing an xz section before (k) and after (l). Orientation axes in (c) and (g) provide scale bars. Nuclei are stained blue (DAPI) and CM are stained green (cTnT). Scale bar: 5 mm (a, d, e) 250 μm, (b, c, f–l).

Figure 3 a: Representative alignment example using a synthetic fiber test image. Input principal axis is (−122°,80°) and input κ is 3.0. Computed local orientations are visualized as small lines, where line color indicates the magnitude in the x (green), y (red), and z (blue) directions. b: Spherical histogram of prescribed fiber orientations. The volume of each ray represents the relative count for each direction and the thick red line shows the principal axis. c: Spherical histogram of the measured local orientations. Scale bars are shown as three-dimensional orientation axes: 200 μm.

We generated porous, large-format (LF) ECTs using 6×10⁶ cells from the same h-iPSC cell/matrix formulation. The LF-ECT geometry was formed using a custom polydimethylsiloxane (PDMS) mold (Fig. 2b). The tissue mold was 21×20.5 mm and 2.5 mm in height. Rectangular PDMS posts, 0.5 mm in width, were evenly distributed over seven columns. In odd columns, two posts, 7 mm in length, were evenly placed. In even columns, a centered 7 mm length post was flanked by two 2.25 mm length posts. The arrangement produced a staggered post pattern, vertically offset by 3.5 mm. After autoclaving, the PDMS mold was coated with 1% Pluronic^® F-127 (Molecular Probes, Eugene, OR, USA) for 1 h and then rinsed with phosphate-buffered saline (PBS). The cell/matrix mixture was manually poured into the mold resulting in gel polymerization, similar to linear ECT generation. LF-ECTs were cultured for 14 or 28 days. The LF-ECT has a diamond mesh shape, where thin tissue “bundles” intersect at “junctions” (Fig. 2c). The current analysis only examined LF-ECT linear bundles. We generated h-iPSC sheet geometry using a similar method, where the PDMS mold had anchors along the perimeter and no internal posts.

The h-iPSC used for the sheet and LF geometries were cultured as described in our previous publication (Masumoto et al., 2016). The h-iPSC used in the linear ECTs, however, were cultured using a Matrigel sandwich protocol, which may produce more mature CM (Zhang et al., 2012).

Cell Transfection for Optogenetic Pacing (OP)

We transfected linear h-iPSC ECTs with an adeno-associated virus at 500 multiplicity of infection (MOI) that delivered the light sensitive ion channel channelrhodopsin-2 (ChIEF), linked to a TdTomato reporter, assisted by Ad-Null virus co-infection (100 MOI). This transfection generated ECTs containing OP enabled CM. Optimal transfection protocols and OP of h-iPSC CM were confirmed in pilot experiments. We then used chronic in vitro OP starting on day 7 using a custom light-emitting diode array powered by an Arduino microcontroller. We recorded intrinsic and maximal OP beat rates on days 7–13 with custom video software and set the chronic pacing rate 0.5 Hz below the maximum OP rate. Chronically paced tissues were compared with control linear h-iPSC ECTs at day 14.

Tissue Preparation and Imaging

We performed immunohistochemistry (IHC) staining and fluorescent confocal microscopy to obtain 3D images of nuclei and CM in whole-mount ECT samples (Figs. 2d–2e). We fixed ECTs for 30 min at room temperature (RT) in 4% paraformaldehyde. Fixed ECTs were washed in 1% Triton X-100/PBS for 1 h at RT and then blocked with 1% Triton X-100/PBS+10% FBS for 1 h at RT. To stain CM, we incubated with primary antibody cardiac troponin T (cTnT, ab45932; Abcam, Cambridge, UK) 1:300 in 1% Triton X-100/PBS+10% FBS+0.2% sodium azide overnight at 4°C followed by a wash with 1% Triton X-100/PBS+10% FBS, incubation with donkey anti-rabbit IgG secondary antibody and then Alexa Fluor 488 conjugated overnight at 4°C. We stained nuclei with 4',6-diamidino-2-phenylindole (DAPI) for 30 min. After IHC staining, we incubated ECTs with 100% glycerol overnight followed by 75% glycerol for 2 h. We then cleared ECTs for 2 h in a solution of 53% benzyl alcohol, 45% glycerol, and 2% DABCO (1,4-diazabicyclo[2,2,2]octane, an antioxidant to preserve dye lifetime) and then changed to fresh clearing solution overnight. This clearing solution has been described previously (Clendenon et al., 2011) and does not require ethanol dehydration. Processed samples were stored in clearing solution at 4°C.

We acquired confocal z-stacks of whole-mount ECT samples with a Nikon ECLIPSE Ti Confocal System (Nikon, Tokyo, Japan) attached to a Nikon Ti-E inverted microscope platform. Images were captured at 1024×1024 pixel density and 3 μm z-step with a 10×0.3 NA objective using Nikon NIS Elements AR software (Nikon). Laser excitation wavelengths were 405 nm for DAPI and 488 nm for cTnT. ECTs were placed in a 10 mm glass bottom dish (Fisher Scientific, Hampton, NH, USA) filled with clearing solution during imaging. Because the refractive indices (RI) of the tissue and imaging medium must be equal for accurate 3D image acquisition, we measured the RI of processed ECTs using optical coherence tomography (Thorlabs, Newton, NJ, USA) and adjusted the clearing solution to match (Tearney et al., 1995). Due to the RI mismatch between the objective (air) and medium (RI=1.52), the actual z-distance scanned inside the tissue is greater than our nominal z-step of 3 μm. To calculate the true distance between slices in the z-stack, we multiplied our nominal z-step by the ratio between the medium and objective RIs (Diaspro et al., 2002). This corrected z-distance (4.56 μm) was used as the voxel z-dimension in our analysis. Both the x- and y-voxel dimensions were 1.24 μm.

Acquisition time for a 100-slice z-stack was ~54 min. We set a maximum acquisition time of 2 h and adjusted the stack size based on this limit. For thinner LF-ECT samples, we were able to acquire the entire tissue. Because linear ECTs were thicker we acquired z-stacks to just beyond the midpoint depth. As our tissues were symmetric, these truncated z-stacks did not affect our measurements. A typical z-stack was 150 slices (679 μm). Images were acquired at 12-bit pixel depth and converted to 8-bit during processing to reduce memory cost. Before applying our orientation method, we corrected image attenuation (Biot et al., 2008) and reduced image noise using a 3D low pass Gaussian blur (SD 3 μm). These preprocessing steps were performed on the cTnT image with ImageJ (NIH, Bethesda, MD, USA).

Linear ECTs and LF-ECT bundles were cylindrical and we computed orientation vectors based on the cylindrical coordinate system defined by the ECT medial axis. To calculate the medial axis, we merged the DAPI and cTnT images and then manually segmented the whole ECT volume with Seg3D (Scientific Computing and Imaging Institute, 2015). We then extracted the medial axis based on the distance transform of the binary image (Schena & Favretto, 2007; Kerschnitzki et al., 2013), smoothed the axis (Garcia, 2010, 2011) and fit a piecewise spline. For a given point x in the ECT, we computed the location p on the medial axis with tangent vector t that satisfied (x−p)·t=0. The axial vector at x was then defined as t, the radial vector as x−p and the circumferential vector as t×(x−p). After collecting the gradient vectors G ^m for a subregion, we transformed to the local ECT coordinate system, where x is the subregion centroid, and then proceeded with the structure tensor analysis.

We implemented our 3D orientation analysis in Matlab (version 8.3.0.532; Mathworks, Natick, MA, USA) using a MacBook Pro laptop (OS×10.9.5; Apple Inc., Cupertino, CA, USA) with a 2.3 GHz Intel Core i7 processor and 16 GB 1600 MHz DDR3 RAM. The raw ND2 files generated during confocal imaging were imported using BioFormats (Linkert et al., 2010). For the 512×512×200 voxel validation images, computation time was 8 min per image, where the gradient operation took 45 s. When generating the kernels for the gradient computation, we used the same voxel spacing as our confocal image data to maintain constant μm units. Results are presented as mean±SD. We performed two-tailed t-tests to examine differences in alignment between ECT groups; p<0.05 was considered significant.

Results

Validation of Principal Axis

We did not find a significant difference between μ measured by our method and the prescribed μ ₀ of the fiber axes for any of the 25 synthetic fiber images (minimum p=0.17). The azimuthal and polar angles measured by our method were within 6° of the prescribed angles (Figs. 3, 4). The mean angle between the measured and prescribed principal axes was 3.2±1.7°. When we compared the x–y–z components of the measured μ and prescribed μ ₀, the two correlated strongly (ρ=0.991, R ²=0.999, Fig. 4), validating that our method accurately measures the principal orientation of 3D images.

Figure 4 Principal axis (μ) validation using synthetic images. a: Actual versus measured azimuthal angle (θ) for each image. b: Actual versus measured polar angles (φ). The actual and measured principal axes were statistically similar for each of test image (n=25 images, p>0.2 for each). c: Actual versus measured principal axes components. The high correlation coefficient (ρ) and low p-value between the actual and measured principal axes indicate strong correlation. The R ² value versus the line of identity shows that the principal axis components are accurately predicted by our algorithm.

Validation of κ

For κ≤4, we did not find any significant differences between the κ measured by our method and the prescribed κ ₀ of the synthetic fiber axes (Fig. 5). However, as κ increased, our algorithm overestimated this metric because we measure many (10-fold) more local orientations than there are fibers. As κ increases, these orientations are often parallel, which generates a greater concentration around μ. Our measurement method may in fact be more accurate of the real CM alignment, as it is intrinsically weighted by fiber length; longer fibers produce more local measurements. Despite differences at high κ, the measured and prescribed values correlated strongly (r=0.996, R ²=0.889, Fig. 5).

Figure 5 Alignment concentration validation. a: Actual versus measured κ for each test image (n=27 images, identity line: gray dash). Significant differences (p<0.05 for individual measures) between the actual and measured κ are indicated in red. Pearson’s coefficient (r) shows strong linear correlation across the complete range of κ values. The R ² coefficient is reduced due to divergence at higher κ. b: The mean κ computed for test images with the same input κ values (κ range from 1 to 10, n=5 images per κ; error bars: SD) were similar (p>0.05 for each κ).

Alternate alignment metrics OOP and α were accurately measured by our method as well (OOP r=0.999, R ²=0.980; α r=0.999, R ²=0.995, Fig. 6). When we compare κ to OOP or α, there is some linear correlation (κ versus OOP r=0.934, κ versus α r=−0.650, Fig. 6). The OOP varies similar to a cosine function of κ, which is expected (equation 5). The cone angle α reduces quickly as κ increases and remains similar for κ>6. These two metrics are therefore less sensitive at high κ values and may obscure differences in tissue organization. For this reason, we use κ as our primary alignment metric.

Figure 6 Alternative alignment metrics show skewed distributions versus κ. a: Actual versus measured orientational order parameter (OOP) for each test image (n=27 images). The OOP ranges from −1 (100% alignment perpendicular to principal axis) to 1 (100% alignment parallel to principal axis). The large Pearson’s coefficient (r) and R ² indicate strong correlation between the actual and measured values. b: Actual versus measured 95% confidence cone angles. The Pearson’s coefficient and R ² value show excellent correlation. c: Correlation between κ and OOP. d: Correlation between κ and the 95% confidence cone angle.

Although our method overestimated κ compared with its expected value, it does result in a robust measurement. When we compared orientations measured for five images designed to have the same κ but different μ, we find that all five do indeed have a common κ (p>0.05), despite it exceeding the expected value (Fig. 5). Therefore, our method produces an alignment metric that is sufficiently precise for comparisons.

CM Alignment in ECTs

Linear and LF-ECTs generated from h-iPSC or embryonic chick cardiac cells underwent comparable gel compaction and began beating in vitro at day 2–3, similar to our previous studies (Tobita et al., 2006; Fujimoto et al., 2011; Masumoto et al., 2016). Both spontaneous beating and beating during optical pacing were synchronized. Preliminary analysis of paraffin sections imaged at 40× showed similar cell density and CM fraction in control versus paced day 14 linear ECTs, suggesting that transfection and chronic pacing do not result in significant cell damage. Whole-mount confocal imaging revealed elliptical cross-sections in both linear ECTs and LF-ECT bundles (Fig. 2). CM within ECTs were predominantly distributed on the outer surface, a pattern observed in our previous experiments, though the roles of metabolic or mechanical stress and biologic regulation establishing this configuration has not been fully elucidated (Tobita et al., 2006; Masumoto et al., 2016). We confirmed this CM distribution with longitudinal paraffin sections and ruled out the possibility of impaired antibody penetration in whole-mount samples.

Whole-mount confocal imaging of cTnT stained ECTs captured CM with sufficient contrast and resolution to observe myocardial organization (Figs. 7–10). The distribution of CM orientations in h-iPSC ECT sheets was significantly less aligned than linear or LF-ECT geometries, but was not entirely isotropic (κ=1.15±0.34, OOP=−0.11±0.07, Figs. 7, 11). At regions close to the sheet edges, CM were more organized due to tension generated by compaction as the sheet pulled away from its anchors. In h-iPSC LF-ECTs, CM predominantly aligned parallel to the bundle long axis (Fig. 8). At day 14, CM were organized around the bundle axis, but the orientations were moderately dispersed (κ=1.62±0.18, OOP=−0.01±0.04, Fig. 8a). Culturing LF-ECTs to day 28 significantly increased CM organization, which resembled linear ECTs derived from embryonic chick cells (κ=3.70±0.56, OOP=0.36±0.08, Fig. 8b). CM in day 14 linear ECTs were strongly aligned to the long axis, for both the h-iPSC and embryonic chick cell sources. Chronic pacing by stimulation of transfected ChIEF did not change CM organization in h-iPSC ECTs (κ=2.97±1.34, OOP=0.23±0.23, nonpaced; κ=3.15±1.17, OOP=0.26±0.21, paced, Figs. 9, 11). Interestingly, alignment in linear h-iPSC ECTs was similar to day 28 h-iPSC LF-ECTs and not significantly different than embryonic chick ECTs (κ=5.16±1.02, OOP=0.53±0.10, Figs. 10, 11). The difference in alignment between day 14 h-iPSC LF and linear ECTs is notable, since both use the same cell source. However, is likely not due to differences in local geometry and may be due to a methods change in our h-iPSC CM induction protocol (Matrigel sandwich; Zhang et al., 2012) used for the h-iPSC linear ECTs. Preliminary functional results show that ECTs generated from h-iPSC differentiated with Matrigel sandwich produce greater active stress and have higher maximum capture rates compared with our previous protocol (Masumoto et al., 2016). The active stress generated by embryonic chick ECTs remains the highest, and is 1.75 times greater than day 14 linear or LF h-iPSC ECT (Clause et al., 2009; Masumoto et al., 2016). Although CM alignment is similar in linear h-iPSC and embryonic chick ECTs, the appearance of CM morphology is qualitatively more mature in the embryonic chick ECTs (Fig. 10). Ongoing research to quantify metrics such as CM size, aspect ratio, or myofiber length, may offer more insight into the differences between h-iPSC and embryonic chick constructs.

Figure 7 a: Local cardiomyocyte (CM) orientations within an human-induced pluripotent stem cells engineered cardiac tissue (ECT) sheet. CM: cardiac troponin T (green). Local CM orientations are visualized as small lines, where line color indicates the magnitude in the x (green), y (red), and z (blue) directions. For clarity, only a subset of orientations is shown. Inset shows region enclosed by dashed lines at 3× scale. b: Spherical histogram of local CM orientations. The volume of each ray represents the relative count for each direction and the thick red line shows the mean CM orientation. κ for this sample was 1.41. Scale bars are shown as three-dimensional orientation axes: 250 μm.

Figure 8 Representative local cardiomyocyte (CM) orientations within human-induced pluripotent stem cells large-format engineered cardiac tissue cultured for (a) 14 days or (b) 28 days. CM: cardiac troponin T (green). Local CM orientations are visualized as small lines, where line color indicates magnitude in the circumferential (green), radial (red), and axial (blue) directions. For clarity, only a subset of the orientations is shown. Insets show regions enclosed by dashed lines at 3× scale. c,d: Spherical histograms display local CM orientations. The volume of each ray represents the relative count for each direction and the thick red line shows the mean CM orientation. κ for the day 14 sample was 1.56 while the day 28 sample was 3.20. Scale bars are shown as three-dimensional orientation axes: 250 μm.

Figure 9 Representative local cardiomyocyte (CM) orientations within (a) day 14 control and (b) day 14 optically paced linear h-iPSC ECTs. CM: cardiac troponin T (green). Local CM orientations are visualized as small lines, where line color indicates magnitude in the circumferential (green), radial (red), and axial (blue) directions. For clarity, only a subset of the orientations is shown. Insets show regions enclosed by dashed lines at 3X scale. c, d: Spherical histograms display local CM orientations. The volume of each ray represents the relative count for each direction and the thick red line shows the mean CM orientation. κ for the control sample was 3.06 while the paced sample was 3.46. Scale bars are shown as 3D orientation axes: 250 μm.

Figure 10 Representative local cardiomyocyte (CM) orientations within a day 14 linear ECT generated from day 10 embryonic chick cardiac cells. CM: cardiac troponin T (green). a: Local CM orientations are displayed as small lines, where line color indicates magnitude in the circumferential (green), radial (red), and axial (blue) directions. For clarity, only a subset of the orientations is shown. Inset shows region enclosed by dashed lines at 3X scale. b: The spherical histogram displays local CM orientations. The volume of each ray represents the relative count for each direction and the thick red line shows the mean CM orientation. κ for this sample was 6.13. Scale bars are shown as 3D orientation axes: 250 μm.

Figure 11 a: Average cardiomyocyte (CM) three-dimensional alignment, κ, computed using our automated gradient method. b: Average CM orientational order parameter (OOP). Tissues were cultured for 14 days unless otherwise noted. Increased CM alignment occurred in large format (LF) and linear engineered cardiac tissue (ECT) geometries versus ECT cell sheet. CM alignment also increased in response to prolonged 28 day culture. Error bars show SD. *p<0.05. iPSC, induced pluripotent stem cells.

CM Local Alignment

We achieve local orientation measurements by dividing an image into subregions. The subregion centroids are evenly distributed over the 3D space, based on the prescribed m _min. For the structure tensor to accurately measure myofiber direction, a subregion should capture a whole cross-section. It is clear that if all subregions are the same size, some will intersect fibers longitudinally, containing only a partial cross-section, and some may be undersized and exist completely inside of a fiber. To account for subregion location and differently sized myofibers, we implemented an adaptive subregion size. We adjust the subregion to minimize the shape parameter derived from the structure tensor eigenvalues. This parameter can discriminate flat, plate-like objects, which occur when a myofiber is longitudinally intersected, and rod-like objects, which are the expected myofiber shape. We found that incorporating this adaptive subregion strategy significantly improved the performance of our method. We repeated our validation test for a subset of κ ₀ values using a fixed subregion size of 10 μm and found that κ measured without the adaptive subregion did not correlate with the prescribed κ ₀. The R ² value for the identity line fit to κ versus κ ₀ was negative (R ²=−0.67), indicating that a horizontal line is a better fit. The adaptive subregion is therefore necessary to achieve accurate orientation and alignment measurements.

DISCUSSION

The majority of studies describing ECT formulations and analytics refer to CM alignment, but most offer only qualitative 2D descriptions (Zimmermann et al., 2000, 2002; Radisic et al., 2004; Huang et al., 2007; Black et al., 2009; Mannhardt et al., 2016). Quantitative results show that geometry, perfusion, culture conditions, and cell population significantly impact CM organization (Tulloch et al., 2011; Bian et al., 2014; Lux et al., 2016; Masumoto et al., 2016). Although variations in cell types, matrix substrates, ECT shape, and alignment metrics make it difficult to compare with our current results, our linear and LF-ECTs follow similar trends of elongation and alignment parallel to the long axis and reduced alignment in h-iPSC-derived tissues compared with late-stage embryonic or neonatal sources (Tulloch et al., 2011; Bian et al., 2014; Hirt et al., 2014; Lux et al., 2016). We previously measured CM alignment in paraffin sections of day 14 linear h-iPSC ECTs using the 2D method described in Karlon et al. (1999). The 3D analysis described here showed less alignment concentration (5.82±1.60 2D versus 2.97±1.34 3D). This difference is due to the limitations of representing 3D tissues in 2D, which neglects geometry and out of plane features (see Masumoto et al., 2016 for the 2D analysis). We noted that optically pacing h-iPSC-derived ECTs above their intrinsic rate for 7 days did not induce a change in CM alignment (Fig. 11). Other studies that used electrical stimulation to accelerate h-iPSC ECT maturation indicated that CM alignment increased (Nunes et al., 2013; Hirt et al., 2014). However, they relied on manual analytic methods conducted on small sample sizes. We only attempted one pacing method and protocol, and alternate methods and protocols may have a greater affect. The direction of the electric field is important in alignment of CM during electrical stimulation (Baumgartner et al., 2015; Maxwell et al., 2016), and a lack of electric field in our optical pacing protocol may therefore diminish the alignment effect. In vivo, chronic pacing effects have been more detrimental than beneficial to maturation. In the embryo, electrical pacing diminished cell proliferation, leading to thin compact myocardium and hypoplasia near the pacing site (Sedmera et al., 1999). In juvenile dogs, prolonged ectopic pacing resulted in CM damage including myofibrillar disarray (Karpawich et al., 1990). Long-term right ventricular apical pacing in adult patients generated electrical dyssynchrony, leading to LV remodeling and dysfunction (Sarkar et al., 2016). Thus, the role of chronic pacing and the underlying cellular and molecular mechanisms affecting maturation remain unknown. Mechanical stretch during tissue culture has been clearly associated with greater alignment and force generation (Tulloch et al., 2011; Godier-Furnemont et al., 2015). We have applied uniaxial stretch to embryonic chick and rat derived ECTs, but have not yet quantified 3D alignment under these conditions (Clause et al., 2009; Ye et al., 2013). During embryonic cardiac development mechanical loading is necessary for differentiation and maturation of the conduction system (Reckova et al., 2003; Sankova et al., 2010), which suggests it may influence electrophysiology of ECTs as well.

Similar to other engineered tissues, CM in our ECTs aligned parallel to the long axis (Costa et al., 2003; Tulloch et al., 2011). Gel compaction may restrict collagen along this axis, which then guides cell alignment. However, earlier studies suggest that cells align collagen, rather than collagen organizes cells. Smooth muscle cells seeded into a prealigned matrix reorganized collagen in response to mechanical loading conditions (Barocas et al., 1998). When exposed to a change in uniaxial loading direction, fibroblasts rapidly remodeled the matrix, even after previous mechanical alignment (Lee et al., 2008). Although the reason for alignment in the unconstrained direction is not fully understood, stress generated during gel compaction is likely a primary factor. The fixed linear ECT ends provide resistance to gel compaction, generating stress in the longitudinal direction. In our LF-ECT, the posts constrain the bundle ends at the junctions, similarly creating longitudinal stress. We do not currently know the magnitude of this stress, but it is sufficient to organize CM. Measuring strains using fiducial markers or modifying the stiffness of our LF-ECT posts could provide quantitative data on the role of mechanical load in the alignment of CM in our ECTs.

Prolonged 28 day culture of h-iPSC LF-ECTs significantly increased CM alignment compared with normal 14 day culture (Fig. 11). The alignment was similar to the linear embryonic chick ECT, which is derived from a more mature embryonic ventricular population. The day 28 tissues were significantly thinner than day 14, indicating continued matrix remodeling throughout this period. Investigations in chick and mouse embryos show that the embryonic myocardial wall is initially isotropic and does not gain a preferential alignment until approximately halfway through development, while cardiac structural and functional maturation continues beyond birth (Tobita et al., 2005; Hirschy et al., 2006). We currently have very limited temporal data on CM organization in the human embryo and fetus. Midgestation (14–27 weeks) fetal hearts analyzed with polarized light microscopy identified organized myofibers as geodesics in a “nested set of warped pretzels” similar to some adult heart models, but did not measure alignment or anisotropy specifically (Jouk et al., 2000). A study using DTMRI to examine alignment at 10, 13, 14, and 19 weeks compared the fetal heart with a 6 day postnatal and adult heart (Mekkaoui et al., 2013). No anisotropic organization was observed until week 19, and it was still significantly less than the day 6 postnatal and adult samples. Our results suggest that the immature h-iPSC CM in our ECTs have a similar, gradual progression of alignment. Although increasing culture time improves tissue organization, it results in a smaller tissue with fewer cells, which may not be sufficient to improve myocardial function after implantation. Therefore, optimal culture requires a balance between tissue organization and cell loss. Accelerating CM maturation with mechanical or electrical stimulation, described above, may overcome this temporal hurdle.

The structural organization of CM within ECTs provided a relevant model to test our 3D alignment method. Our method is based on a generic analysis of 3D histologic images, and is therefore applicable to investigations beyond cardiac tissue engineering. Tissue clearing and confocal microscopy have been used to reconstruct whole heart tissue architecture in avian embryos (Miller et al., 2005; Mao et al., 2013) and the adult guinea pig and mouse (Smith et al., 2008). Combining these techniques with our method could quantify 3D cardiac structure and investigate the development of myocardial organization in the embryo. In neurobiology, tracking alignment of growing neurites can further our understanding of neuronal development and improve the efficacy of neuroregenerative therapies (Mahoney et al., 2005; Portera-Cailliau et al., 2005; Kofron et al., 2010). Clinically, quantifying neurite orientations in neurodegenerative diseases such as Alzheimer’s could establish criteria to differentiate from normal tissue (Saxena & Caroni, 2007). Studies on arterial mechanics measure smooth muscle, elastin, and collagen organization to develop constitutive models of normal and diseased states, which suggest that structural orientation may have a role in growth and rupture (Elbischger et al., 2004; Eriksson et al., 2009; Wilson et al., 2012). Disorders in skeletal muscle associated with myofiber disarray, including Duchenne muscular dystrophy, a degenerative disease, and McArdle disease, a metabolic myopathy, require analysis of local orientations to assess muscle damage (Wang et al., 2015; Krag et al., 2016). Investigators have adapted 2D alignment techniques to examine these areas, but the complex tissue structures are better represented in 3D (Elbischger et al., 2004; Smith et al., 2008; Mitchel et al., 2013; Gunther et al., 2015). Currently, directed imaging modalities or image segmentation methods are applied for 3D analysis (Gasser et al., 2012; Wang et al., 2015; Singh et al., 2016). Our method presents an automated 3D approach that accurately computes local orientations and alignment at the cellular level.

Limitations

Our method determines orientation based on image gradient and is therefore dependent on contrast at object boundaries. The gradient computation is also sensitive to image noise, which we reduce by applying a low pass Gaussian filter. The imaging z-step is a critical factor as it determines vertical resolution and image acquisition time. A small z-step requires more slices while a large z-step does not capture vertical orientations. We performed tests using our synthetic fiber images, generated using different distances between slices, and found that as values increased above 6 μm, the measured μ became skewed toward the vertical axis. However, z-steps below 3 μm required image acquisition times of more than 3 h per ECT. Extensive imaging time can result in photobleaching or evaporation of imaging medium. Therefore, we chose 3 μm as our z-step during imaging, which corresponded to a z-distance of 4.56 μm after correcting for RI mismatch between the objective and medium. To obtain 3D images of ECTs, we used a clearing solution of benzoyl alcohol and glycerol, which preserves tissue structure better than benzyl alcohol/benzyl benzoate solutions, but may still cause distortion (Clendenon et al., 2011). We selected the benzoyl alcohol and glycerol clearing method due to its simplicity and because it does not require preceding dehydration. As demonstrated in a recent comparative study, different clearing protocols produce varying degrees of tissue transparency and fluorescent preservation, and the appropriate method should be selected based on the research objective (Kolesova et al., 2016). Despite clearing, attenuation occurred during imaging and limited our scan depth. We applied a prospective attenuation correction to correct this issue (Biot et al., 2008). Alternatively, hardware-based depth correction during imaging, such as increasing excitation intensity and photodetector gain as a function of z-position, could avoid attenuation effects (Smith et al., 2008). Finally, the computation time for our method increases linearly with the number of voxels, which may become impractical for large images. Converting to a compiled code such as C or C++ or using parallel computing would significantly reduce computation time.

CONCLUSIONS

CM alignment in ECTs is associated with faster conduction (Bian et al., 2014), more rapid Ca²⁺ cycling (Feinberg et al., 2012; Khan et al., 2015), and increased stress production (Black et al., 2009). These enhanced functions are desirable for implantable constructs used in MI repair. The method presented and validated in this study provides an accurate measurement of CM organization in 3D. The use of tissue clearing and confocal microscopy enables whole-tissue imaging at a resolution sufficient to calculate CM orientations, while incorporating an adaptive subregion size accounts for different scales and produces greater accuracy. Alignment metrics such as κ and OOP are robust and detect differences in tissue organization. A faster confocal system and generating compiled code would make the time cost for larger tissue imaging and analysis more reasonable. Although we tested our method with in vitro cardiac tissues, it is applicable to a variety of research areas requiring quantitative 3D structure.