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., Reference Banglmaier, Sander and Vandevord2015), ocular tissue (Avila & Bueno, Reference Avila and Bueno2015), neurite growth (Singh et al., Reference Singh, Negi, Laezza, Papadakis and Labate2016), vasculature (Gasser et al., Reference Gasser, Ogden and Holzapfel2006), and the heart (Streeter et al., Reference Streeter, Spotnitz, Patel, Ross and Sonnenblick1969; Gilbert et al., Reference Gilbert, Benson, Li and Holden2007). 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., Reference Matsuda, Takahashi, Nariai, Ito, Takatani, Fujio and Azuma2005; Nguyen et al., Reference Nguyen, Tinney, Yuan, Roussel, El-Baz, Giridharan, Keller and Sethu2013) and endothelial cells (EC) elongate and align in the direction of steady flow or cyclic stretch (Galbraith et al., Reference Galbraith, Skalak and Chien1998; Kaunas et al., Reference Kaunas, Nguyen, Usami and Chien2005). Genetic errors associated with abnormal myofiber alignment and fibrosis form the basis of a range of clinically important cardiomyopathies (Carrier et al., Reference Carrier, Bonne, Bahrend, Yu, Richard, Niel, Hainque, Cruaud, Gary, Labeit, Bouhour, Dubourg, Desnos, Hagege, Trent, Komajda, Fiszman and Schwartz1997; Akdis et al., Reference Akdis, Brunckhorst, Duru and Saguner2016; Finsterer & Stollberger, Reference Finsterer and Stollberger2016) and are associated with abnormal myocardial structure and function in congenital heart disease (Hahn et al., Reference Hahn, Lauriol, Thul, Behnke-Hall, Logeswaran, Schanzer, Bogurcu, Garvalov, Zenker, Gelb, Von Gerlach, Kandolf, Kontaridis and Schranz2015).
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., Reference Gilbert, Benson, Li and Holden2007). Establishing this myoarchitecture during normal cardiac morphogenesis requires complex temporal and spatial myocardial maturation (Sedmera et al., Reference Sedmera, Pexieder, Vuillemin, Thompson and Anderson2000; Tobita et al., Reference Tobita, Garrison, Liu, Tinney and Keller2005). 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., Reference Waldman, Nosan, Villarreal and Covell1988; Omens et al., Reference Omens, May and Mcculloch1991; Costa et al., Reference Costa, Takayama, Mcculloch and Covell1999; Ashikaga et al., Reference Ashikaga, Van Der Spoel, Coppola and Omens2009). 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, Reference Panfilov and Keener1993; Punske et al., Reference Punske, Taccardi, Steadman, Ershler, England, Valencik, Mcdonald and Litwin2005; Hooks et al., Reference Hooks, Trew, Caldwell, Sands, Legrice and Smaill2007; Dhein et al., Reference Dhein, Seidel, Salameh, Jozwiak, Hagen, Kostelka, Hindricks and Mohr2014). Structural organization is established during early embryonic development, where myofiber orientations proceed through several iterations (Tobita et al., Reference Tobita, Garrison, Liu, Tinney and Keller2005). When the embryonic heart is exposed to abnormal mechanical loading, myofiber organization is disrupted, indicating biomechanical regulation (Sedmera et al., Reference Sedmera, Pexieder, Vuillemin, Thompson and Anderson2000; Tobita et al., Reference Tobita, Garrison, Liu, Tinney and Keller2005) and producing altered developmental trajectories (Burggren & Reyna, Reference Burggren and Reyna2011). 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., Reference Wickline, Verdonk, Wong, Shepard and Miller1992; Sosnovik et al., Reference Sosnovik, Wang, Dai, Wang, Aikawa, Novikov, Rosenzweig, Gilbert and Wedeen2009; Rutherford et al., Reference Rutherford, Trew, Sands, Legrice and Smaill2012). In the hypertrophic heart, the transmural helical variation of CM is reduced and tissue organization is diminished (Roberts & Ferrans, Reference Roberts and Ferrans1975; Schmitt et al., Reference Schmitt, Fedarava, Falkenberg, Rothaus, Bodhey, Reischauer, Kozerke, Schnackenburg, Westermann, Lunkenheimer, Anderson, Berger and Kuehne2009; Hill et al., Reference Hill, Simon, Valdez-Jasso, Zhang, Champion and Sacks2014). Myocardial fibrosis in patients with hypertrophic cardiomyopathy is associated with increased electrical heterogeneity and, importantly, significant risk of sudden, arrhythmic death (Kadish et al., Reference Kadish, Shinnar, Moore, Levine, Balke and Spear1988; Valderrabano et al., Reference Valderrabano, Lee, Ohara, Lai, Fishbein, Lin, Karagueuzian and Chen2001; Ripplinger et al., Reference Ripplinger, Li, Hadley, Chen, Rothenberg, Lombardi, Wickline, Marian and Efimov2007; Muthappan & Calkins, Reference Muthappan and Calkins2008).
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., Reference Mozaffarian, Benjamin, Go, Arnett, Blaha, Cushman, Das, De Ferranti, Despres, Fullerton, Howard, Huffman, Isasi, Jimenez, Judd, Kissela, Lichtman, Lisabeth, Liu, Mackey, Magid, Mcguire, Mohler, Moy, Muntner, Mussolino, Nasir, Neumar, Nichol, Palaniappan, Pandey, Reeves, Rodriguez, Rosamond, Sorlie, Stein, Towfighi, Turan, Virani, Woo, Yeh and Turner2016). 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., Reference Bolli, Chugh, D’amario, Loughran, Stoddard, Ikram, Beache, Wagner, Leri, Hosoda, Sanada, Elmore, Goichberg, Cappetta, Solankhi, Fahsah, Rokosh, Slaughter, Kajstura and Anversa2011; Sanganalmath & Bolli, Reference Sanganalmath and Bolli2013; Masumoto et al., Reference Masumoto, Ikuno, Takeda, Fukushima, Marui, Katayama, Shimizu, Ikeda, Okano, Sakata and Yamashita2014, Reference Masumoto, Nakane, Tinney, Yuan, Ye, Kowalski, Minakata, Sakata, Yamashita and Keller2016; Fisher et al., Reference Fisher, Doree, Mathur and Martin-Rendon2015). 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., Reference Ribeiro, Ang, Fu, Rivas, Mohamed, Higgs, Srivastava and Pruitt2015; Li et al., Reference Li, Mkrtschjan, Lin and Russell2016), micropatterned coatings (Feinberg et al., Reference Feinberg, Alford, Jin, Ripplinger, Werdich, Sheehy, Grosberg and Parker2012; Tsang et al., Reference Tsang, Annabi, Ercole, Zhou, Karst, Li, Haynes, Evans, Thissen, Khademhosseini and Forsythe2015), and tissue geometry (Bian et al., Reference Bian, Jackman and Bursac2014; Nguyen et al., Reference Nguyen, Badie, Mcspadden, Pedrotty and Bursac2014) all influence CM size, shape, and alignment. In addition, mechanical and electrical conditioning may accelerate tissue organization (Nunes et al., Reference Nunes, Miklas, Liu, Aschar-Sobbi, Xiao, Zhang, Jiang, Masse, Gagliardi, Hsieh, Thavandiran, Laflamme, Nanthakumar, Gross, Backx, Keller and Radisic2013; Stoppel et al., Reference Stoppel, Kaplan and Black2016). 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., Reference Chaudhuri, Nguyen, Rangayyan, Walsh and Frank1987; Petroll et al., Reference Petroll, Cavanagh, Barry, Andrews and Jester1993). Fitting an ellipse model to cell boundaries can quantify individual orientations, but is highly dependent on the segmentation method (Bray et al., Reference Bray, Adams, Geisse, Feinberg, Sheehy and Parker2010; Xu et al., Reference Xu, Beyazoglu, Hefner, Gurkan and Demirci2011; Nectow et al., Reference Nectow, Kilmer and Kaplan2014). In cases where ellipse models are not appropriate, segmented images have been processed using the Hough transform (Bayan et al., Reference Bayan, Levitt, Miller, Kaplan and Georgakoudi2009) or eigenanalysis of the 2nd-order moment tensor (Vader et al., Reference Vader, Kabla, Weitz and Mahadevan2009) 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., Reference Chaudhuri, Kundu and Sarkar1993; Hong et al., Reference Hong, Wan and Jain1998; Karlon et al., Reference Karlon, Hsu, Li, Chien, Mcculloch and Omens1999; Grosberg et al., Reference Grosberg, Alford, Mccain and Parker2011; Sun et al., Reference Sun, Bloom and Zaman2015). Our group has applied this method to quantify myofiber alignment in the embryonic heart and in ECTs (Tobita et al., Reference Tobita, Garrison, Liu, Tinney and Keller2005; Masumoto et al., Reference Masumoto, Nakane, Tinney, Yuan, Ye, Kowalski, Minakata, Sakata, Yamashita and Keller2016). 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., Reference Hsu, Muzikant, Matulevicius, Penland and Henriquez1998) and polarized light microscopy (Tower & Tranquillo, Reference Tower and Tranquillo2001) 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., Reference Wu, Rajwa, Filmer, Hoffmann, Yuan, Chiang, Sturgis and Robinson2003) 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 Reference Bredfeldt, Liu, Conklin, Keely, Mackie and Eliceiria , Reference Bredfeldt, Liu, Pehlke, Conklin, Szulczewski, Inman, Keely, Nowak, Mackie and Eliceiri2014 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., Reference Frangi, Niessen, Vincken and Viergever1998; Daniels et al., Reference Daniels, Ter Haar Romeny, Rubbens and Van Assen2006). This method is improved when applying an adaptive filter size, though it does significantly increase computation time (Daniels et al., Reference Daniels, Ter Haar Romeny, Rubbens and Van Assen2006). Other methods rotate local subregions to maximize the z-axis 2nd-order moment (Reuze et al., Reference Reuze, Coatrieux, Luo and Dillenseger1993), 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., Reference Robb, Wirjadi and Schladitz2007; Wirjadi et al., Reference Wirjadi, Schladitz, Rack and Breuel2009), or compute the distances to the fiber boundary in all 26 directions to identify the major inertia axis (Altendorf & Jeulin, Reference Altendorf and Jeulin2009; Altendorf et al., Reference Altendorf, Decenciere, Jeulin, De Sa Peixoto, Deniset-Besseau, Angelini, Mosser and Schanne-Klein2012). 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, Reference Weickert1998; Puspoki et al., Reference Puspoki, Storath, Sage and Unser2016). 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., Reference Zhao, Butters, Zhang, Pullan, Legrice, Sands and Smaill2012 a, Reference Zhao, Krueger, Seemann, Meng, Zhang, Dossel, Legrice and Smaill2012 b) and in high-resolution 3D cardiovascular magnetic resonance images of explanted rat hearts (Gilbert et al., Reference Gilbert, Benoist, Benson, White, Tanner, Holden, Dobrzynski, Bernus and Radjenovic2012; Bernus et al., Reference Bernus, Radjenovic, Trew, Legrice, Sands, Magee, Smaill and Gilbert2015). 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, Reference Tabor and Rokita2007; Kersh et al., Reference Kersh, Zysset, Pahr, Wolfram, Larsson and Pandy2013; Larsson et al., Reference Larsson, Luisier, Kersh, Dall’ara, Zysset, Pandy and Pahr2014). 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, Reference Tabor and Rokita2007; Kersh et al., Reference Kersh, Zysset, Pahr, Wolfram, Larsson and Pandy2013; Larsson et al., Reference Larsson, Luisier, Kersh, Dall’ara, Zysset, Pandy and Pahr2014).
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. (Reference Karlon, Hsu, Li, Chien, Mcculloch and Omens1999), 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.
$$G_{x} {\equals}h_{x}\, {\ast}\,I,\,G_{y} {\equals}h_{y} \,{\ast}\,I,\,G_{z} {\equals}h_{z} \,{\ast}\,I,$$
$$\!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 1, … , 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 (|λ 1|≤|λ 2|≤|λ 3|) and corresponding eigenvectors (e 1, e 2, e 3) of T. The eigenvector e 1, 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., Reference Fisher, Lewis and Embleton1987). 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., Reference Fisher, Lewis and Embleton1987).
$${\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., Reference Fisher, Lewis and Embleton1987). 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.
$$\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 1, … , 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, Reference Fisher1995). The most common measurement in 2D, however, is the orientational order parameter (OOP, equation 5; Drew et al., Reference Drew, Eagleson, Baldo, Parker and Grosberg2015). 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 (ψ).
$${\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., Reference Fisher, Lewis and Embleton1987). 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 (λ 3) 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., Reference Fisher, Lewis and Embleton1987).
$${{\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. (Reference Fisher, Lewis and Embleton1987, §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., Reference Chen, Wei, Newstadt, Degraef, Simmons and Hero2015). We varied the principal axis (μ 0) or concentration parameter (κ 0) 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 μ 0 and κ 0. To test μ, κ 0 was fixed to 3.0 and we varied the azimuthal (θ) and polar (φ) angles of μ 0 between −180° and 180° and 10° and 90°, respectively, and generated one image for each μ 0 (25 total images). To test κ, μ 0 was fixed to (135°, 45°) and we varied κ 0 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 0). 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., Reference Fisher, Lewis and Embleton1987). We then examined correlation between the set of measured principal axes (M=[μ 1, … , μ 25]) and prescribed principal axes (M 0=[μ 0,1, … , μ 0,25]), where a coefficient (ρ) of 1 indicates perfect correlation (Fisher et al., Reference Fisher, Lewis and Embleton1987). We also computed the R 2 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 0 had equal concentrations (Fisher, Reference Fisher1986). 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 κ 0 values. In addition, we computed the R 2 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 κ 0 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 μ 0 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 κ 0, we tested if κ for all of the measured orientation distributions were the same (Fisher, Reference Fisher1986). 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., Reference Clause, Tinney, Liu, Keller and Tobita2006), embryonic rat (Fujimoto et al., Reference Fujimoto, Clause, Liu, Tinney, Verma, Wagner, Keller and Tobita2011), and human-induced pluripotent stem cells (h-iPSC)-derived cardiac cells (Masumoto et al., Reference Masumoto, Nakane, Tinney, Yuan, Ye, Kowalski, Minakata, Sakata, Yamashita and Keller2016) 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., Reference Masumoto, Nakane, Tinney, Yuan, Ye, Kowalski, Minakata, Sakata, Yamashita and Keller2016). Embryonic chick CM were harvested from day 10 white Leghorn chick embryos. Collected cells (~3.0×106 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 NaHCO3, 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 CO2 incubator (37°C, 5% CO2) 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−5 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×106 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., Reference Masumoto, Nakane, Tinney, Yuan, Ye, Kowalski, Minakata, Sakata, Yamashita and Keller2016). 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., Reference Zhang, Klos, Wilson, Herman, Lian, Raval, Barron, Hou, Soerens, Yu, Palecek, Lyons, Thomson, Herron, Jalife and Kamp2012).
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., Reference Clendenon, Young, Ferkowicz, Phillips and Dunn2011) 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., Reference Tearney, Brezinski, Southern, Bouma, Hee and Fujimoto1995). 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., Reference Diaspro, Federici and Robello2002). 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., Reference Biot, Crowell, Hofte, Maurin, Vernhettes and Andrey2008) 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, Reference Schena and Favretto2007; Kerschnitzki et al., Reference Kerschnitzki, Kollmannsberger, Burghammer, Duda, Weinkamer, Wagermaier and Fratzl2013), smoothed the axis (Garcia, Reference Garcia2010, Reference Garcia2011) 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., Reference Linkert, Rueden, Allan, Burel, Moore, Patterson, Loranger, Moore, Neves, Macdonald, Tarkowska, Sticco, Hill, Rossner, Eliceiri and Swedlow2010). 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 μ 0 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 μ 0, the two correlated strongly (ρ=0.991, R 2=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 2 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 κ 0 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 2=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 2 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 2=0.980; α r=0.999, R 2=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 2 indicate strong correlation between the actual and measured values. b: Actual versus measured 95% confidence cone angles. The Pearson’s coefficient and R 2 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., Reference Tobita, Liu, Janczewski, Tinney, Nonemaker, Augustine, Stolz, Shroff and Keller2006; Fujimoto et al., Reference Fujimoto, Clause, Liu, Tinney, Verma, Wagner, Keller and Tobita2011; Masumoto et al., Reference Masumoto, Nakane, Tinney, Yuan, Ye, Kowalski, Minakata, Sakata, Yamashita and Keller2016). 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., Reference Tobita, Liu, Janczewski, Tinney, Nonemaker, Augustine, Stolz, Shroff and Keller2006; Masumoto et al., Reference Masumoto, Nakane, Tinney, Yuan, Ye, Kowalski, Minakata, Sakata, Yamashita and Keller2016). 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., Reference Zhang, Klos, Wilson, Herman, Lian, Raval, Barron, Hou, Soerens, Yu, Palecek, Lyons, Thomson, Herron, Jalife and Kamp2012) 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., Reference Masumoto, Nakane, Tinney, Yuan, Ye, Kowalski, Minakata, Sakata, Yamashita and Keller2016). 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., Reference Clause, Tinney, Liu, Keller and Tobita2009; Masumoto et al., Reference Masumoto, Nakane, Tinney, Yuan, Ye, Kowalski, Minakata, Sakata, Yamashita and Keller2016). 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 κ 0 values using a fixed subregion size of 10 μm and found that κ measured without the adaptive subregion did not correlate with the prescribed κ 0. The R 2 value for the identity line fit to κ versus κ 0 was negative (R 2=−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., Reference Zimmermann, Fink, Kralisch, Remmers, Weil and Eschenhagen2000, Reference Zimmermann, Schneiderbanger, Schubert, Didie, Munzel, Heubach, Kostin, Neuhuber and Eschenhagen2002; Radisic et al., Reference Radisic, Park, Shing, Consi, Schoen, Langer, Freed and Vunjak-Novakovic2004; Huang et al., Reference Huang, Khait and Birla2007; Black et al., Reference Black, Meyers, Weinbaum, Shvelidze and Tranquillo2009; Mannhardt et al., Reference Mannhardt, Breckwoldt, Letuffe-Breniere, Schaaf, Schulz, Neuber, Benzin, Werner, Eder, Schulze, Klampe, Christ, Hirt, Huebner, Moretti, Eschenhagen and Hansen2016). Quantitative results show that geometry, perfusion, culture conditions, and cell population significantly impact CM organization (Tulloch et al., Reference Tulloch, Muskheli, Razumova, Korte, Regnier, Hauch, Pabon, Reinecke and Murry2011; Bian et al., Reference Bian, Jackman and Bursac2014; Lux et al., Reference Lux, Andree, Horvath, Nosko, Manikowski, Hilfiker-Kleiner, Haverich and Hilfiker2016; Masumoto et al., Reference Masumoto, Nakane, Tinney, Yuan, Ye, Kowalski, Minakata, Sakata, Yamashita and Keller2016). 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., Reference Tulloch, Muskheli, Razumova, Korte, Regnier, Hauch, Pabon, Reinecke and Murry2011; Bian et al., Reference Bian, Jackman and Bursac2014; Hirt et al., Reference Hirt, Boeddinghaus, Mitchell, Schaaf, Bornchen, Muller, Schulz, Hubner, Stenzig, Stoehr, Neuber, Eder, Luther, Hansen and Eschenhagen2014; Lux et al., Reference Lux, Andree, Horvath, Nosko, Manikowski, Hilfiker-Kleiner, Haverich and Hilfiker2016). We previously measured CM alignment in paraffin sections of day 14 linear h-iPSC ECTs using the 2D method described in Karlon et al. (Reference Karlon, Hsu, Li, Chien, Mcculloch and Omens1999). 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., Reference Masumoto, Nakane, Tinney, Yuan, Ye, Kowalski, Minakata, Sakata, Yamashita and Keller2016 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., Reference Nunes, Miklas, Liu, Aschar-Sobbi, Xiao, Zhang, Jiang, Masse, Gagliardi, Hsieh, Thavandiran, Laflamme, Nanthakumar, Gross, Backx, Keller and Radisic2013; Hirt et al., Reference Hirt, Boeddinghaus, Mitchell, Schaaf, Bornchen, Muller, Schulz, Hubner, Stenzig, Stoehr, Neuber, Eder, Luther, Hansen and Eschenhagen2014). 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., Reference Baumgartner, Halbach, Krausgrill, Maass, Srinivasan, Sahito, Peinkofer, Nguemo, Muller-Ehmsen and Hescheler2015; Maxwell et al., Reference Maxwell, Wagner and Davis2016), 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., Reference Sedmera, Grobety, Reymond, Baehler, Kucera and Kappenberger1999). In juvenile dogs, prolonged ectopic pacing resulted in CM damage including myofibrillar disarray (Karpawich et al., Reference Karpawich, Justice, Cavitt and Chang1990). Long-term right ventricular apical pacing in adult patients generated electrical dyssynchrony, leading to LV remodeling and dysfunction (Sarkar et al., Reference Sarkar, Tilkar, Jain, Mondal, Sarkar and Modi2016). 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., Reference Tulloch, Muskheli, Razumova, Korte, Regnier, Hauch, Pabon, Reinecke and Murry2011; Godier-Furnemont et al., Reference Godier-Furnemont, Tiburcy, Wagner, Dewenter, Lammle, El-Armouche, Lehnart, Vunjak-Novakovic and Zimmermann2015). 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., Reference Clause, Tinney, Liu, Keller and Tobita2009; Ye et al., Reference Ye, Yuan, Li, Cooper, Tinney and Keller2013). During embryonic cardiac development mechanical loading is necessary for differentiation and maturation of the conduction system (Reckova et al., Reference Reckova, Rosengarten, Dealmeida, Stanley, Wessels, Gourdie, Thompson and Sedmera2003; Sankova et al., Reference Sankova, Machalek and Sedmera2010), 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., Reference Costa, Lee and Holmes2003; Tulloch et al., Reference Tulloch, Muskheli, Razumova, Korte, Regnier, Hauch, Pabon, Reinecke and Murry2011). 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., Reference Barocas, Girton and Tranquillo1998). When exposed to a change in uniaxial loading direction, fibroblasts rapidly remodeled the matrix, even after previous mechanical alignment (Lee et al., Reference Lee, Holmes and Costa2008). 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., Reference Tobita, Garrison, Liu, Tinney and Keller2005; Hirschy et al., Reference Hirschy, Schatzmann, Ehler and Perriard2006). 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., Reference Jouk, Usson, Michalowicz and Grossi2000). 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., Reference Mekkaoui, Porayette, Jackowski, Kostis, Dai, Sanders and Sosnovik2013). 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., Reference Miller, Thompson, Bigelow, Gittinger, Trusk and Sedmera2005; Mao et al., Reference Mao, Gribble, Pertsov and Shi2013) and the adult guinea pig and mouse (Smith et al., Reference Smith, Matiukas, Zemlin and Pertsov2008). 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., Reference Mahoney, Chen, Tan and Saltzman2005; Portera-Cailliau et al., Reference Portera-Cailliau, Weimer, De Paola, Caroni and Svoboda2005; Kofron et al., Reference Kofron, Liu, Lopez-Fagundo, Mitchel and Hoffman-Kim2010). Clinically, quantifying neurite orientations in neurodegenerative diseases such as Alzheimer’s could establish criteria to differentiate from normal tissue (Saxena & Caroni, Reference Saxena and Caroni2007). 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., Reference Elbischger, Bischof, Regitnig and Holzapfel2004; Eriksson et al., Reference Eriksson, Kroon and Holzapfel2009; Wilson et al., Reference Wilson, Baek and Humphrey2012). 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., Reference Wang, Zhang, Wasala, Duan and Yao2015; Krag et al., Reference Krag, Pinos, Nielsen, Brull, Andreu and Vissing2016). Investigators have adapted 2D alignment techniques to examine these areas, but the complex tissue structures are better represented in 3D (Elbischger et al., Reference Elbischger, Bischof, Regitnig and Holzapfel2004; Smith et al., Reference Smith, Matiukas, Zemlin and Pertsov2008; Mitchel et al., Reference Mitchel, Martin and Hoffman-Kim2013; Gunther et al., Reference Gunther, Gunther, Schneiders, Rupp and Blesch2015). Currently, directed imaging modalities or image segmentation methods are applied for 3D analysis (Gasser et al., Reference Gasser, Gallinetti, Xing, Forsell, Swedenborg and Roy2012; Wang et al., Reference Wang, Zhang, Wasala, Duan and Yao2015; Singh et al., Reference Singh, Negi, Laezza, Papadakis and Labate2016). 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., Reference Clendenon, Young, Ferkowicz, Phillips and Dunn2011). 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., Reference Kolesova, Capek, Radochova, Janacek and Sedmera2016). Despite clearing, attenuation occurred during imaging and limited our scan depth. We applied a prospective attenuation correction to correct this issue (Biot et al., Reference Biot, Crowell, Hofte, Maurin, Vernhettes and Andrey2008). 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., Reference Smith, Matiukas, Zemlin and Pertsov2008). 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., Reference Bian, Jackman and Bursac2014), more rapid Ca2+ cycling (Feinberg et al., Reference Feinberg, Alford, Jin, Ripplinger, Werdich, Sheehy, Grosberg and Parker2012; Khan et al., Reference Khan, Xu, Hua, Johnson, Belevych, Janssen, Gyorke, Guan and Angelos2015), and increased stress production (Black et al., Reference Black, Meyers, Weinbaum, Shvelidze and Tranquillo2009). 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.
Acknowledgments
This work was supported by the Kosair Charities Pediatric Heart Research Endowment (to B.B.K.). The authors also thank Dr. Patricia Soucy and the staff of the Speed School of Engineering confocal imaging core for imaging assistance.
