In a properly contracting cardiac muscle, many different subcellular structures are organized into an intricate architecture. While it has been observed that this organization is altered in pathological conditions, the relationship between length-scales and architecture has not been properly explored. In this work, we utilize a variety of architecture metrics to quantify organization and consistency of single structures over multiple scales, from subcellular to tissue scale as well as correlation of organization of multiple structures. Specifically, as the best way to characterize cardiac tissues, we chose the orientational and co-orientational order parameters (COOPs). Similarly, neonatal rat ventricular myocytes were selected for their consistent architectural behavior. The engineered cells and tissues were stained for four architectural structures: actin, tubulin, sarcomeric z-lines, and nuclei. We applied the orientational metrics to cardiac cells of various shapes, isotropic cardiac tissues, and anisotropic globally aligned tissues. With these novel tools, we discovered: (1) the relationship between cellular shape and consistency of self-assembly; (2) the length-scales at which unguided tissues self-organize; and (3) the correlation or lack thereof between organization of actin fibrils, sarcomeric z-lines, tubulin fibrils, and nuclei. All of these together elucidate some of the current mysteries in the relationship between force production and architecture, while raising more questions about the effect of guidance cues on self-assembly function. These types of metrics are the future of quantitative tissue engineering in cardiovascular biomechanics.

## Introduction

The mechanical pumping of the heart is indisputably a function of the multiscale architecture of the myocardium [13]. In cardiac tissue engineering, recapitulating this architecture over multiple length-scales has been one of the barriers for the utilization of stem-cell derived cardiomyocytes for fixing the heart and for the creation of in vitro platforms for testing cardiotoxicity of pharmaceuticals [4,5]. For both of these applications, it is increasingly evident that the architecture of the scaffold affects maturity of cells, gene expression levels, and various functional outputs [612]. One affected function is the cells' ability to self-assemble into a tissue with optimized contractile properties. Learning to leverage the cells' self-assembly capability would allow for tissues to be engineered with minimal input. However, to accomplish this it is necessary to quantify the tissue architecture over many scales and identify the essential cues that promote recapitulation of native tissue structure in engineered cardiac tissue.

During normal development, the heart macrostructure evolves from a tube into a consistent anatomy composed of four chambers wrapped with helical fibers [79,1316]. Similarly, normal primary single cardiomyocytes interrogated for a response to geometrical extracellular matrix (ECM) cues will produce qualitatively consistent structures [17,18]. However, it is unclear if this consistency holds overall length-scales, and no one has been able to reliably quantify this. Previously, quantifying orientation of architecture in cells has been done using various computational algorithms [1926]. Also, prior studies used a variety of metrics to analyze their results that were based on different statistical distributions, such as Gaussian [12,20,21,24,26,27]. Yet, organizations of structures in cardiac cells are not often normally distributed and one of the best suited parameters for this task is or based on the orientational order parameter (OOP). While initially developed for the field of liquid crystals, these have now been extensively used to quantify organization in biology [2833]. Still, there is no standard method to measure the architecture of various structures (nuclei, actin fibrils, microtubules, etc.) over multiple scales to ensure that in vitro tissues recapitulate the condition of interest. Quantifying this architecture would also allow us to answer the important question of how local and global organizations differ in their effect on contractile function. To properly build in vitro replicas, it is also essential to quantify the correlation between organizations of different subcellular structures. This is evident from observing the orientation of α-actinin structures (protein in the z-lines of the myofibril) and the actin fibrils. In development, these two evolve from parallel, punctuate, in stress fibers to perpendicular z-lines in myofibrils [34]. In sum, to control the mechanics of engineered heart tissues, it is essential to understand multiscale consistency, architecture, and the correlation of various structures in cardiac muscle.

The heart architecture has been extensively studied for decades [1,12,16,33,3537]. However, the complex nature of the structures that form the heart makes it challenging to quantitatively describe architecture in a consistent manner. This consistency is essential as the field begins to rely more on cardiac cells and tissues engineered from stem cells [5,38]. Currently, these are often assembled as two-dimensional (2D) monolayers, which can be thought of as an essential building block of the heart because the myocardium is organized with sheets of myofibrils arranged in parallel [39]. In order to both understand and rebuild the three-dimensional (3D) myocardium, it is necessary to fully understand the 2D sheets. To this end, we demonstrate how a collection of order parameters [28,29,33,37] can be used to elucidate cardiac cell and tissue architecture over multiple scales that span from subcellular to tissue scale. The qualitative observation of consistent self-assembly as a function of cellular shape was tested over multiple subcellular length-scales and was found to have implications on these cells contractile function. Further, we quantify the patches of organization spontaneously assembled in engineered isotropic tissues and show the correlation between orientation of actin fibrils and z-lines, microtubules, and nuclei. In all, this work illustrates the challenge of recapitulating the nonsimplistic multiscale architecture of cardiac tissues, but also its importance in creating engineered tissues with proper mechanical function.

## Methods

Standard methods were used to prepare cells and tissues for structure characterization [17,37]. Here, we describe them briefly.

### Substrate Preparation.

Glass coverslips (Fisher Scientific Company, Hanover Park, IL) were sonicated in a 95% ethanol solution for 30 min. These coverslips were then spin-coated in polydimethylsiloxane mixed at a 10:1 ratio with curing agent (PDMS, Ellsworth Adhesives, Germantown, WI). The coverslips with PDMS coating were cured at 65 °C overnight.

### Extracellular Matrix Printing.

Stamp designs were drawn using adobe illustrator software (Adobe Systems, Inc., San Jose, CA). The designs were etched into 5 in. × 5 in. chrome with soda-lime glass masks by a third-party vendor (FrontRange Photo Mask Co., Palmer Lake, CO). Silicon wafers were made through SU-8 deposition using the glass masks in the Bio-Organic Nanofabrication Facility (University of California, Irvine, Irvine, CA). In order to create the stamps, silicon wafers were covered in 60–80 g of PDMS. It was then cured at 65 °C overnight and peeled from the wafers. Thereafter, the square patterned regions were cut out and stored for use as stamps.

The stamps were sonicated and incubated for 1 h with 0.1 mg/mL drops of fibronectin (Fisher Scientific Company, Hanover Park, IL) and placed on top of the patterned faces. Prior to applying fibronectin to the coverslips, the coverslips were ultraviolet-ozone-treated. Then, the coverslips were incubated in 1% pluronic acid solution (5 g pluronic F-127, Sigma Aldrich, Inc., Saint Louis, MO, dissolved in 500 mL sterile water) for 10 min in order to block cell adhesion between regions of fibronectin. Coverslips designated for isotropic samples were placed with the PDMS side facing downward on 300 μL drops of 0.05 mg/mL of fibronectin for 10 min. After these steps, coverslips were rinsed in phosphate-buffered saline (PBS) (Life Technologies, Carlsbad, CA) three times and kept in PBS until cell seeding.

### Neonatal Rat Ventricular Myocytes (NRVM) Harvest and Seeding.

All the animal protocols were reviewed and approved by the University of California, Irvine, Institutional Animal Care and Use Committee (Protocol No. 2013-3093). Neonatal Sprague-Dawley rats (Charles River Laboratories, Wilmington, MA), 1–3 days postpartum, were sprayed with 95% ethanol and decapitated. The hearts were rapidly removed and trimmed in Hank's balanced salt solution (HBSS) (Life Technologies, Carlsbad, CA). Once all the hearts (ten per harvest) were dissected, they were incubated at 4 °C overnight (12 h) in a 1 mg/mL trypsin solution (Sigma-Aldrich, Inc., Saint Louis, MO) dissolved in HBSS. The trypsin solution was then removed, and the tissue was neutralized in warmed M199 culture medium (Invitrogen, Carlsbad, CA) supplemented with 10% heat-inactivated fetal bovine serum, 10 mM HEPES, 20 mM glucose, 2 mM L-glutamine (Life Technologies, Carlsbad, CA), 1.5 μM vitamin B-12, and 50 U/ml penicillin (Sigma-Aldrich, Inc., Saint Louis, MO). The media were removed, and the tissue was dissociated through several washes of 1 mg/mL collagenase dissolved in HBSS. The collagenase cell solutions were then centrifuged at 1200 rpm for 10 min. The supernatant was aspirated, and the cells were resuspended in chilled HBSS. The HBSS cell solution was centrifuged at 1200 rpm for 10 min. The supernatant was aspirated and the cells were resuspended in warm 10% FBS M199. The cells were purified in three preplate steps, including 45 min incubations in two different cell culture flasks and a 40 min incubation in a third cell culture flask (BD Biosciences, San Diego, CA) within a tissue culture incubator. Cells were counted using a disposable hemocytometer (Fisher Scientific, Waltham, MA) and seeded at a density of 1 × 106 cells per 3 mL for anisotropic and isotropic coverslips. Cells were seeded at a density of 1.5 × 105 or 2.5 × 105 cells per 3 mL for single cell shape coverslips.

The samples for every experiment were cultured in standard incubator conditions. The dead cells were washed off the samples 24 h after seeding when media were refreshed with 10% FBS M199 media. The media were replaced by 2% FBS M-199 media 48 h after seeding to maintain the cardiac myocytes without rapidly increasing the fibroblast population.

### Fixing and Staining.

After 72 h in culture, the samples were fixed in 4% paraformaldehyde (Fisher Scientific Company, Hanover Park, IL) supplemented with 0.001% Triton X-100 (Sigma Aldrich, Inc., Saint Louis, MO). The samples were each stained for a combination of four or three of the following: actin (Alexa Fluor 488 Phalloidin, Life Technologies, Carlsbad, CA), sarcomeric α-actinin (mouse monoclonal anti-α-actinin, Sigma Aldrich, Inc., Saint Louis, MO), nuclei (4′,6′-diaminodino-2-phenylinodole, Life Technologies, Carlsbad, CA), fibronectin (polyclonal rabbit antihuman fibronectin, Sigma Aldrich, Inc., Saint Louis, MO), and alpha-tubulin (mouse monoclonal anti-α-tubulin or chicken polyclonal anti-α-tubulin; Abcam plc., San Francisco, CA). Secondary staining was done using the appropriate pairs including: tetramethylrhodamine-conjugated goat antimouse IgG antibody (Alexa Fluor 633 Goat Anti-Mouse or Alexa Fluor 750 Goat Anti-Mouse, Life Technologies, Carlsbad, CA), tetramethylrhodamine-conjugated goat antichicken IgG antibody (Alexa-Fluor 488 Goat Anti-Chicken, Life Technologies, Carlsbad, CA), and goat antirabbit IgG antibodies (Alexa Fluor 750 Goat Antirabbit, Life Technologies, Carlsbad, CA). Coverslips with stained tissues were mounted onto glass cover slides (VWR, Radnor, PA) with ProLong Gold Antifade Reagent (Life Technologies, Carlsbad, CA) to prevent stain bleaching. Nail polish (Electron Microscopy Sciences, Hatfield, PA) was then used as sealant along the edge of each coverslip. The sealant was allowed to dry overnight.

### Imaging.

Cover slides with immunostained samples were imaged on an IX-83 inverted motorized microscope (Olympus America, Center Valley, PA) mounted with a digital CCD camera ORCA-R2 C10600-10B (Hamamatsu Photonics, Shizuoka Prefecture, Japan) using an UPLFLN 40 × oil immersion objective (Olympus America, Center Valley, PA). Ten fields of view were acquired for each anisotropic and isotropic sample.

### Image Analysis: Tissues.

Image processing was done using ImageJ software (NIH, Bethesda, MD) and customized matlab software (MathWorks, Natick, MA) as previously described [12,33,37] to determine orientation of structures and metrics. One of the metrics, the orientation order parameter has been commonly used for over 40 yr [28,40]. More recently, it has been utilized to describe organization for biological applications [2933]. The parameter can be derived with a variety of methods including calculating coefficients of the Fourier series and by taking the eigenvalues of the structure tensor [17,41]. The orientational order parameter is commonly labeled as OOP [33,41], f2D [31], or S [28]. By convention, the order parameter ranges from 0 to 1, which in the structure tensor method can be accomplished by taking the deviatoric portion. This method is briefly summarized below for an image that contained a subcellular structure with orientations $p(x,y)$
$T=〈2[pi,xpi,xpi,xpi,ypi,xpi,ypi,ypi,y]−[1001]〉$
(1)
where the operator 〈 〉 is the average over all vectors $p$ (i.e., over parameter i in Eq. (1)). Then, orientational order parameter (OOP) is simply the maximum eigenvalue of the average tensor T and ranges from 0 to 1
$OOP=max(eigenvalue(〈T〉))$
(2)
The co-orientational order parameter (COOP) for two structures with orientations $p(x,y)$ and $q(x,y)$ is also simply the maximum eigenvalue of the appropriate tensor which is calculated as follows:
$fi,x=pi,xqi,x+pi,yqi,y and fi,y=pi,xqi,x−pi,xqi,x$
(3)

$TPQ=〈2[fi,xfi,xfi,xfi,yfi,xfi,yfi,yfi,y]−[1001]〉$
(4)

$COOP=max(eigenvalue(TPQ))$
(5)
Also, the uncorrelated COOP and correlated COOP were calculated
(6)

(7)

If the COOP is identical to the correlated COOP, then the two structures are as correlated as possible. Conversely, if the COOP is identical to the uncorrelated COOP, the two structures are as uncorrelated as possible [37]. The COOP ranges from 0 to 1, and its limits (i.e., uncorrelated and correlated COOPs) can also vary from 0 to 1 with the corresponding orientational order in the two structures.

Each field of view consisted of an area of 3.5 × 104 μm2 at 40× magnification. Global organization was defined by ten fields of view at 40× magnification, totaling an area of approximately 0.35 mm2, while the total area of the tissue is on the order of 1 cm2.

### Image Analysis: Single Cells.

Image analysis of single cells was performed using previously created customized matlab codes. These codes were used to determine subcellular structure orientation and calculate the COOP [37]. Grid squares were placed over the entire image, and the director, i.e., average orientation, was calculated in each grid square. The area of each cell was approximately 1250 μm2, which is equivalent to a length-scale of 35 μm. Therefore, the possible sizes of the grid squares used to analyze single cardiomyocytes could range from 1 μm to 35 μm. We chose the large length-scale (21 μm) by determining the largest scale that contained a minimum of two orientation angles within two grid squares for all the cell shapes. Pairwise comparisons were used to quantify the consistency of cardiomyocytes of same shape and aspect ratio (AR).

### Nuclei Detect and Voronoi Code.

Nuclei detection was performed using a custom matlab code written to detect and segment fluorescently stained cell nuclei from the background of an image. An intensity threshold was applied to the image in order to generate a binary mask representing each nucleus as a body of white pixels and the background as black. The boundary of each nucleus was then calculated as a set of points in 2D space, allowing automated, quantitative characterization of nuclear size, shape, and orientation [4244].

Following automated detection and segmentation, each cell nucleus within an image is assigned a centroid, calculated from its best-fit ellipse [44]. These centroids are then used to generate a Voronoi diagram, which partitions the image into spaces consisting of all the points closest to a specific nuclear centroid [45]. Each Voronoi approximates the boundary of the area influenced by each nucleus. Then, the single nucleus orientation value is expanded into a vector field spanning its entire Voronoi space and compared to the actin or tubulin orientational vector field falling within the space.

### Statistics.

To determine statistical significance, one-way analysis of variance was used with the Holm-Sidak post-hoc test, which is commonly used for pairwise comparison of experimental groups. Significance was considered for an unadjusted p-value less than the critical level, which accounts for the number of comparisons.

## Results

Characterizing cardiac muscle architecture is inherently a multiscale problem. To address this, the consistency of self-assembly of single cardiomyocytes over multiple length-scales was considered first.

### Consistency of Single Cell Assembly.

There are two questions of interest when considering single cells: (1) how does consistency of self-assembly change in a multiscale way and (2) how does extracellular matrix (ECM) cues affect consistency of self-assembly? In cell shapes with parallel fibers (i.e., rectangles with high aspect ratio), it is possible to evaluate consistency based on the variance of a simple order parameter. For cells that have unaligned myofibrils (i.e., squares), the orientational order parameter (OOP) would be always low, and the variance would be independent of consistency. To overcome this, the co-orientational order parameter (COOP) was utilized to evaluate consistency of self-assembly in cells of various shapes. This provided a method to evaluate consistency as a function of scale in single cardiomyocytes by culturing NRVMs on fibronectin islands of rectangles, ovals, and triangles of various aspect ratios (Table 1; Fig. 1(a)). The COOP was calculated to quantify consistency of actin fibrils and sarcomeric z-lines for a variety of subcellular length-scales (Fig. 1(b)). For all the shapes, the COOP varied from low values at small length-scales to high values at large length-scales (Fig. 2). To summarize, the consistency data of actin fibrils and sarcomeric z-lines were plotted for all the shapes at a small length-scale of ∼1 μm (Fig. 2(a); and Fig. S1A, which is available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection) and at a large length-scale of ∼21 μm (Fig. 2(b); and Fig. S1B, which is available under the “Supplemental Materials” tab). Consistency of self-assembly at the small scale (∼1 μm) was highly dependent on the aspect ratio of the ECM island (Fig. 2(a); and Fig. S1A, which is available under the “Supplemental Materials” tab). Furthermore, consistency of actin fibrils was significantly different between circles (ovals with aspect ratio of 1) and various shapes with an aspect ratio of ∼1 (Fig. 2(a)). Ovals with aspect ratios of 13.49 and 15.06 were significantly different than rectangles with an aspect ratio of 14.20 (Fig. 2(a)). Also, at this small scale, the consistency of actin fibrils settles at an aspect ratio of ∼6, which correlates with the known relationship between force and aspect ratio [46]. The average COOP at a scale of ∼1 μm for shapes with aspect ratios (ARs) >5.5 was 0.87 and 0.39 for actin (Fig. 2(a)) and sarcomeric z-lines (Fig. S1A is available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection), respectively. Except for the circle, the shaped NRVMs were consistent at length-scales of ∼21 μm (area = 441μm2) and above, with an average COOP of 0.98 and 0.87 for actin (Fig. 2(b)) and sarcomeric z-lines (Fig. S1B is available under the “Supplemental Materials” tab), respectively. The oval with an aspect ratio of 1, i.e., a circle, was shown to be less consistent than other shapes at large scales (Fig. 2(b); and Fig. S1B, which is available under the “Supplemental Materials” tab) because it has the added challenge of shape rotational symmetry, but it has also been inconsistent qualitatively [17]. For all the other shapes, the qualitative observation of consistent cell self-assembly of actin fibrils was confirmed by this data for large length-scales. This data illustrated that it is possible to elucidate the self-assembly consistency in single cells, but tissue architecture needed to be characterized over a larger range of length-scale.

### Multiscale Order in Cardiac Tissues.

To address this question, two types of cardiac tissues, anisotropic globally aligned (Figs. 3(a)3(d)-(i)) and isotropic (Figs. 3(a)3(d)-(ii))) tissues made from primary cardiomyocytes (NRVMs), were used. These tissues were stained for nuclei (Fig. 3(a)), sarcomeric z-lines (Fig. 3(b)), actin fibrils (Fig. 3(c)), and microtubules (Fig. 3(d)). The global orientational order was previously measured by calculating the OOP for ten fields of view (A = 3.5 mm2) [12,33], but this only provided the orientational order at one length-scale. We were able to measure the orientational order at smaller length-scales by calculating the OOP for a smaller number of fields of view, and then averaging the results for field of view sets within the tissue sample. Further, each field of view image was divided into multiple quadrants, and the OOP was calculated and averaged for each of the smaller image sets. With this method, it was possible to calculate OOP for areas over 4 orders of magnitude (Fig. 3(e)) that spans from the subcellular to tissue length-scales. Unsurprisingly, for structures that were aligned extremely well globally, such as actin, tubulin, and nuclei, the OOP remains high at smaller subcellular length-scales. It was interesting that the z-lines, which are in general not as well aligned as actin fibrils, maintained the exact same orientational order at smaller subcellular area scales (Fig. 3(e)). This implied that the difference between order of the actin fibrils and sarcomeric z-lines was a local phenomenon, which makes sense since these two structures are part of the same structure—the sarcomere.

As shown previously [12,33], isotropic tissues were found to have a very low global, A = 3.5 × 105 μm2, OOP (Fig. 3(f)). Comparatively, at smaller subcellular length-scales, the OOP increased significantly. For A < 2 × 103 μm2, the tested quadrant was similar in size to or smaller than a single cell, which implied that mostly each quadrant would have a single nuclei. Thus, mathematically an OOP ≈ 1 was expected for the nuclei at these areas. This was confirmed with the nuclei OOP steadily increasing until at the small subcellular scale it was reasonably close to one (Fig. 3(f)). The same logic could not be applied to the other structures considered as they fully populate the cell area. However, a similar increasing trend was observed for actin, tubulin, and z-lines. Actin fibrils and microtubules were significantly better organized at the small subcellular area scales, although neither reached the maximum demonstrated in globally aligned tissues. Conversely, the sarcomeric z-line organization reached the maximum demonstrated in globally aligned tissues (dashed lowest line compared to the triangular data point at A ≈ 1.5 × 103 μm2 in Fig. 3(f)). Collectively, these results proved the qualitative observation that in isotropic tissues, there is a spontaneous formation of well-organized patches (solid lines—Fig. 3(f)). The size of these spontaneously organized patches was approximated by fitting the Hill function to the isotropic OOP data and extracting the EC50 from the fit (Fig. 4(g)). The average characteristic area of the organized patches (based on actin, tubulin, and z-lines) was 1.4 × 104 μm2. Still, the different behaviors of z-lines compared to actin and tubulin at the various length-scales raise the question of how the structure orientations are correlated.

### Correlation of Subcellular Structures in Cardiac Tissues.

To address this, the COOP was utilized to compare orientations of subcellular structure pairs within the same tissue. As a demonstration, the actin fibrils were paired with sarcomeric z-lines (Fig. 4(a)), nuclei (Fig. 4(b)-(i)), and microtubules (Fig. 4(b)-(i)–(ii)). When comparing two structures, the COOP has a range that can be achieved which changes as a function of the OOP of the two structures [37]. As a result, it was not possible to resolve the correlation for globally aligned tissues where in the best case 0.44 ± 0.10 < COOP (aligned Actin, z-lines) < 0.50 ± 0.12. However, in isotropic tissues the range was wide enough to produce statistical significance (Fig. 4(c)). The lower limit was the uncorrelated COOP (circles—Fig. 4(c))—at this limit, the two structures had no detectable correlation in orientations. Conversely, the upper limit was the correlated COOP (diamonds—Fig. 4(c))—at this limit, the two structures were maximally correlated. The actual parameter value (squares—Fig. 4(c)) was then tested for statistical significance against both limits. Unsurprisingly, as they are part of the same sarcomeric structure, actin fibrils were found to be maximally correlated to the z-lines (red—Fig. 4(c)). Microtubules, which are part of the same cytoskeleton, were found to be somewhat correlated to actin fibril orientation (yellow—Fig. 4(c)). Another interesting comparison was the nuclei orientation and the localized cytoskeleton. To compare the nuclei orientation to cytoskeleton fibrils, it was necessary to define some area of cytoskeleton relevant to each nucleus. In this case, the tissue was broken-up by using Voronoi spaces (Fig. 4(d)). While these were not the cell boundaries, they were a reasonable first estimate of the area of influence for each nucleus. To compare a single vector of the nucleus orientation to the many actin or tubulin fiber vectors in the same space. Each spatial position that had an actin vector was also assigned a nucleus vector value based on its Voronoi space (Figs. 4(e) and 4(f)). The resultant COOP values showed that the nucleus orientation was maximally uncorrelated with respect to the actin fibril orientation (Fig. 4(c-(iii))). This method was used to compare tubulin and nuclei orientations, which were also significantly uncorrelated (Fig. 4(c)-(iv)). When interpreting this data, it was important to also consider that the nucleus eccentricity was found to be 0.82 ± 0.01 for isotropic tissues and 0.87 ± 0.01 for globally aligned, which was a statistically significant difference with N = 6 tissues.

## Discussion

In this work, we illustrated the use of two order parameters in characterizing multiscale organization of engineered cardiac tissues from the subcellular to tissue scale: The classical orientational order parameter [28,29,33] and the co-orientational order parameter [37]. Indeed, there were other metrics that could have been used to summarize orientation properties from the data obtained through image analysis that are based on a variety of functions. For example, there are different statistical metrics based on different distributions that were used to quantify orientation data [20,21,26,27] or multiple image correlation methods to compare images [25,4750]. However, the OOP and the COOP were chosen because they are mathematically appropriate to interpret data from biological tissues. Also, the OOP and COOP can be used on all types of distributions and, by now, have been widely characterized for use in engineered cardiac tissues [4,9,12,33,36,37]. However, it is important to note that the orientational order parameter is not capable of differentiating between isotropic distributions and ones with multiple principal orientations. By its similar mathematical construction, the COOP also cannot be used to describe multiple correlation types within a single tissue [37]. Thus, if the goal is to capture properties of highly detailed, multipeaked distributions, the COOP and OOP are not the appropriate parameters. These methods and software could have also been used on immunostained images of fibronectin as well [51], however, examining ECM was beyond the scope of this paper.

There have been many models of single cell self-assembly, but they have mostly relied on a qualitative view of cytoskeleton architecture for validation [17,5254]. The ability to quantify consistency displayed in this work provides a tool for future models not only to be validated quantitatively but also to be assessed at multiple subcellular length-scales (Fig. 1(b)). Indeed, by considering architecture over multiple length-scales, we discovered that the highest actin fibril consistency at smaller subcellular length-scales was observed with high aspect ratios, e.g., elongated cells, without sharp corners of the ECM island. Therefore, at small length-scales in order for actin fibrils to consistently self-assemble, it is more advantageous for cells to be a smooth elongated shape. The consistency settled at a very high value for the small subcellular scale of 1 μm when the aspect ratio is above ∼6 (Fig. 2(a)). This has an interesting correlation with the traction force microscopy data showing a maximum force produced by cells with an aspect ratio of 6–7 [46]. Thus, in order to produce reliably high forces the cell needs to be assembled consistently within areas that range 1400 μm2 (Fig. 2), or conversely, given proper guidance, cardiomyocytes can be assembled very consistently with areas of 3 orders of magnitude in scale.

The emergent force production in tissues also highly depends on the multiscale architecture of cardiac muscle [12,36]. It has been shown that for engineered locally anisotropic parquet tiles of ∼6 × 104μm2, with a global isotropic order, the force production is easily predicted based on a net force vector addition model [36]. Conversely, the force production in engineered isotropic tissues (Fig. 3(b)) is significantly lower than expected [12,33,36]. In this work, we quantitatively confirmed that engineered isotropic tissues spontaneously assemble organize patches of 1.4 × 104 μm2, which is the same order of magnitude as tissues composed of parquet tiles proven to be similar to globally aligned tissues (Fig. 3(a)) in their local force generation [36]. It will be the subject of future studies to find out if this is the result of a substantial difference between local organization of 1.4 × 104 μm2 and 6 × 104μm2, or more likely, the result of providing engineered tissue guidance, which influences the cascade of organization. In this respect, understanding the correlation of orientations provides a powerful tool in elucidating the changes in these tissues at the local subcellular scale.

The correlation of sarcomeric components, actin fibrils, and z-lines, can be used to represent the development state [34]. In the primary engineered tissues tested in this work, the correlation is maximally tight (Fig. 4(c)-(i)), but for immature cardiomyocytes with stress fibers like those made from stem-cell derived cardiomyocytes [4], the correlation would be expected to be significantly lower. The absence of correlation between the nucleus and the cytoskeleton components, i.e., actin (Fig. 4(c)-(iii)) and tubulin (Fig. 4(c)-(iv)), was very surprising because in an aligned tissue the nucleus is thought to be oriented due to the forces applied to it by the other structures of the cytoskeleton [55]. There are a few possible explanations for this unexpected result. First, the nuclei in the engineered isotropic tissue were on average less elliptical, which means their directions were less accurately determined. Second, the area of the influence of each nucleus was based on a mathematical estimate, but it might be more applicable to define it based on a biological relevance. For example, excluding the area around the nucleus, which might be wrapped in a net of fibrils, might provide a different result. Nonetheless, the COOP provides a powerful tool to investigate these issues and to quantify correlations of orientations even when it requires a careful selection of structures for comparison in order to provide relevance to mechanics.

## Conclusions

The quantitative analysis of primary engineered cardiac cells and tissues architecture illuminated some important properties: (1) The self-assembly and force development mechanisms are optimized at the same minimal aspect ratio and (2) there are tight correlations in orientations between different subcellular structures within a cytoskeleton, which is important to force generation. The analysis also exposed new mysteries. Indeed, engineered isotropic tissues spontaneously assemble organized patches that are comparable in size to engineered tiles, which deepens the mystery of why the isotropic tissues do not produce more force [12,36]. This might be explained by the loose correlation of nucleus orientation in the cytoskeleton, which in turn might explain the gene expression differences between isotropic and aligned tissues. Both of these are subjects for future investigation. Looking forward, it is clear that such architectural metrics are essential for elucidating the mechanical behavior of engineered cardiac tissues.

## Acknowledgment

We thank Danny Baldo, Jr., Zachery R. Robinson, Tuyetnhi Le, Marlen Tagle Rodriguez, and Priyanka Ganesh for participating in gathering single cell data used in this work. We are indebted to Professor Susan Rafelski and Dr. Irina Mueller for the suggestion to look at microtubules in our tissues. This work was partially funded by National Science Foundation EAGER Award No. 1338609.

A.G. designed research; N.K.D., N.E.J., and J.Q.C. performed experiments; N.K.D, N.E.J., J.Q.C, and A.G. performed data analysis; and N.K.D, N.E.J., J.Q.C, and A.G. wrote the article.

## References

References
1.
Helm
,
P.
,
Beg
,
M. F.
,
Miller
,
M. I.
, and
Winslow
,
R. L.
,
2005
, “
Measuring and Mapping Cardiac Fiber and Laminar Architecture Using Diffusion Tensor MR Imaging
,”
,
1047
(
1
), pp.
296
307
.
2.
Helm
,
P. A.
,
Younes
,
L.
,
Beg
,
M. F.
,
Ennis
,
D. B.
,
Leclercq
,
C.
,
Faris
,
O. P.
,
McVeigh
,
E.
,
Kass
,
D.
,
Miller
,
M. I.
, and
Winslow
,
R. L.
,
2006
, “
Evidence of Structural Remodeling in the Dyssynchronous Failing Heart
,”
Circ. Res.
,
98
(
1
), pp.
125
132
.
3.
Vendelin
,
M.
,
Bovendeerd
,
P.
,
Engelbrecht
,
J.
, and
Arts
,
T.
,
2002
, “
Optimizing Ventricular Fibers: Uniform Strain or Stress, But Not ATP Consumption, Leads to High Efficiency
,”
Am. J. Physiol.
,
283
(3), pp.
1072
1081
.
4.
Sheehy
,
S. P.
,
Pasqualini
,
F.
,
Grosberg
,
A.
,
Park
,
S. J.
,
Aratyn-Schaus
,
Y.
, and
Parker
,
K. K.
,
2014
, “
Quality Metrics for Stem Cell-Derived Cardiac Myocytes
,”
Stem Cell Rep.
,
2
(
3
), pp.
282
294
.
5.
Feric
,
N. T.
, and
,
M.
,
2016
, “
Maturing Human Pluripotent Stem Cell-Derived Cardiomyocytes in Human Engineered Cardiac Tissues
,”
,
96
, pp.
110
134
.
6.
Cohn
,
J. N.
,
Ferrari
,
R.
, and
Sharpe
,
N.
,
2000
, “
Cardiac Remodeling—Concepts and Clinical Implications: A Consensus Paper From an International Forum on Cardiac Remodeling
,”
J. Am. Coll. Cardiol.
,
35
(
3
), pp.
569
582
.
7.
Bursac
,
N.
,
Parker
,
K. K.
,
Iravanian
,
S.
, and
Tung
,
L.
,
2002
, “
Cardiomyocyte Cultures With Controlled Macroscopic Anisotropy—A Model for Functional Electrophysiological Studies of Cardiac Muscle
,”
Circ. Res.
,
91
(
12
), pp.
E45
E54
.
8.
Bursac
,
N.
,
Loo
,
Y. H.
,
Leong
,
K.
, and
Tung
,
L.
,
2007
, “
Novel Anisotropic Engineered Cardiac Tissues: Studies of Electrical Propagation
,”
Biochem. Biophys. Res. Commun.
,
361
(
4
), pp.
847
853
.
9.
McCain
,
M. L.
,
Sheehy
,
S. P.
,
Grosberg
,
A.
,
Goss
,
J. A.
, and
Parker
,
K. K.
,
2013
, “
Recapitulating Maladaptive, Multiscale Remodeling of Failing Myocardium on a Chip
,”
,
110
(
24
), pp.
9770
9775
.
10.
Spach
,
M. S.
,
Heidlage
,
J. F.
,
Dolber
,
P. C.
, and
Barr
,
R. C.
,
2000
, “
Electrophysiological Effects of Remodeling Cardiac Gap Junctions and Cell Size: Experimental and Model Studies of Normal Cardiac Growth
,”
Circ. Res.
,
86
(
3
), pp.
302
311
.
11.
Spach
,
M. S.
,
Heidlage
,
J. F.
,
Barr
,
R. C.
, and
Dolber
,
P. C.
,
2004
, “
Cell Size and Communication: Role in Structural and Electrical Development and Remodeling of the Heart
,”
Heart Rhythm
,
1
(
4
), pp.
500
515
.
12.
Feinberg
,
A. W.
,
Alford
,
P. W.
,
Jin
,
H.
,
Ripplinger
,
C. M.
,
Werdich
,
A. A.
,
Sheehy
,
S. P.
,
Grosberg
,
A.
, and
Parker
,
K. K.
,
2012
, “
Controlling the Contractile Strength of Engineered Cardiac Muscle by Hierarchal Tissue Architecture
,”
Biomaterials
,
33
(
23
), pp.
5732
5741
.
13.
Larsen
,
W. J.
,
1997
, “
Development of the Heart
,”
Human Embryology
,
W. B. Saunders
,
New York
, pp.
131
155
.
14.
Forouhar
,
A. S.
,
Liebling
,
M.
,
Hickerson
,
A.
,
,
A.
,
Tsai
,
H.-J.
,
Hove
,
J. R.
,
Fraser
,
S. E.
,
Dickinson
,
M. E.
, and
Gharib
,
M.
,
2006
, “
The Embryonic Vertebrate Heart Tube is a Dynamic Suction Pump
,”
Science
,
312
(
5774
), pp.
751
753
.
15.
Liebling
,
M.
,
Forouhar
,
A. S.
,
Wolleschensky
,
R.
,
Zimmermann
,
B.
,
Ankerhold
,
R.
,
Fraser
,
S. E.
,
Gharib
,
M.
, and
Dickinson
,
M. E.
,
2006
, “
Rapid Three-Dimensional Imaging and Analysis of the Beating Embryonic Heart Reveals Functional Changes During Development
,”
Dev. Dyn.
,
235
(
11
), pp.
2940
2948
.
16.
Streeter
,
D. D.
, Jr.
,
Spotnitz
,
H. M.
,
Patel
,
D. P.
,
Ross
,
J.
, and
Sonnenblick
,
E. H.
,
1969
, “
Fiber Orientation in the Canine Left Ventricle During Diastole and Systole
,”
Circ. Res.
,
24
(
3
), pp.
339
347
.
17.
Grosberg
,
A.
,
Kuo
,
P. L.
,
Guo
,
C. L.
,
Geisse
,
N. A.
,
Bray
,
M. A.
,
,
W. J.
,
Sheehy
,
S. P.
, and
Parker
,
K. K.
,
2011
, “
Self-Organization of Muscle Cell Structure and Function
,”
PLoS Comput. Biol.
,
7
(
2
), p.
e1001088
.
18.
Sheehy
,
S. P.
,
Grosberg
,
A.
, and
Parker
,
K. K.
,
2012
, “
The Contribution of Cellular Mechanotransduction to Cardiomyocyte Form and Function
,”
Biomech. Model. Mechanobiol.
,
11
(
8
), pp.
1227
1239
.
19.
Bray
,
M.-A. P.
,
,
W. J.
,
Geisse
,
N. A.
,
Feinberg
,
A. W.
,
Sheehy
,
S. P.
, and
Parker
,
K. K.
,
2010
, “
Nuclear Morphology and Deformation in Engineered Cardiac Myocytes and Tissues
,”
Biomaterials
,
31
(
19
), pp.
5143
5150
.
20.
Karlon
,
W. J.
,
Hsu
,
P. P.
,
Li
,
S.
,
Chien
,
S.
,
McCulloch
,
A. D.
, and
Omens
,
J. H.
,
1999
, “
Measurement of Orientation and Distribution of Cellular Alignment and Cytoskeletal Organization
,”
Ann. Biomed. Eng.
,
27
(
6
), pp.
712
720
.
21.
Dunn
,
G. A.
, and
Brown
,
A. F.
,
1986
, “
Alignment of Fibroblasts on Grooved Surfaces Described by a Simple Geometric Transformation
,”
J. Cell Sci.
,
83
, pp.
313
340
.
22.
Wax
,
A.
,
Yang
,
C. H.
,
Backman
,
V.
,
,
K.
,
Boone
,
C. W.
,
Dasari
,
R. R.
, and
Feld
,
M. S.
,
2002
, “
Cellular Organization and Substructure Measured Using Angle-Resolved Low-Coherence Interferometry
,”
Biophys. J.
,
82
(
4
), pp.
2256
2264
.
23.
Molitoris
,
J. M.
,
Paliwal
,
S.
,
Sekar
,
R. B.
,
Blake
,
R.
,
Park
,
J.
,
Trayanova
,
N. A.
,
Tung
,
L.
, and
Levchenko
,
A.
,
2016
, “
Precisely Parameterized Experimental and Computational Models of Tissue Organization
,”
Integr. Biol.
,
8
(
2
), pp.
230
242
.
24.
Kartasalo
,
K.
,
Pölönen
,
R.-P.
,
Ojala
,
M.
,
,
J.
,
Lekkala
,
J.
,
Aalto-Setälä
,
K.
, and
Kallio
,
P.
,
2015
, “
CytoSpectre: A Tool for Spectral Analysis of Oriented Structures on Cellular and Subcellular Levels
,”
BMC Bioinf.
,
16
(
1
), pp.
1
23
.
25.
Zaritsky
,
A.
,
Obolski
,
U.
,
,
Z.
,
Gan
,
Z.
,
Du
,
Y.
,
Schmid
,
S. L.
, and
Danuser
,
G.
,
2016
, “
Decoupling Global and Local Relations in Biased Biological Systems
,”
bioRxiv
, epub.
26.
Morrill
,
E. E.
,
Tulepbergenov
,
A. N.
,
Stender
,
C. J.
,
Lamichhane
,
R.
,
Brown
,
R. J.
, and
Lujan
,
T. J.
,
2016
, “
A Validated Software Application to Measure Fiber Organization in Soft Tissue
,”
Biomech. Model. Mechanobiol.
, epub.
27.
Davidson
,
P.
,
Bigerelle
,
M.
,
Bounichane
,
B.
,
Giazzon
,
M.
, and
Anselme
,
K.
,
2010
, “
Definition of a Simple Statistical Parameter for the Quantification of Orientation in Two Dimensions: Application to Cells on Grooves of Nanometric Depths
,”
Acta Biomater.
,
6
(
7
), pp.
2590
2598
.
28.
Hamley
,
I. W.
,
2007
,
Introduction to Soft Matter: Synthetic and Biological Self-Assembling Materials
,
Wiley
,
Hoboken, NJ
.
29.
Volfson
,
D.
,
Cookson
,
S.
,
Hasty
,
J.
, and
Tsimring
,
L. S.
,
2008
, “
Biomechanical Ordering of Dense Cell Populations
,”
,
105
(
40
), pp.
15346
15351
.
30.
Sun
,
J. G.
,
Tang
,
J.
, and
Ding
,
J. D.
,
2009
, “
Cell Orientation on a Stripe-Micropatterned Surface
,”
Chin. Sci. Bull.
,
54
(
18
), pp.
3154
3159
.
31.
Umeno
,
A.
, and
Ueno
,
S.
,
2003
, “
Quantitative Analysis of Adherent Cell Orientation Influenced by Strong Magnetic Fields
,”
IEEE Trans. Nanobiosci.
,
2
(
1
), pp.
26
28
.
32.
Balachandran
,
K.
,
Alford
,
P. W.
,
Wylie-Sears
,
J.
,
Goss
,
J. A.
,
Grosberg
,
A.
,
Bischoff
,
J.
,
Aikawa
,
E.
,
Levine
,
R. A.
, and
Parker
,
K. K.
,
2011
, “
Cyclic Strain Induces Dual-Mode Endothelial-Mesenchymal Transformation of the Cardiac Valve
,”
,
108
(
50
), pp.
19943
19948
.
33.
Grosberg
,
A.
,
Alford
,
P. W.
,
McCain
,
M. L.
, and
Parker
,
K. K.
,
2011
, “
Ensembles of Engineered Cardiac Tissues for Physiological and Pharmacological Study: Heart on a Chip
,”
Lab Chip
,
11
(
24
), pp.
4165
4173
.
34.
Du
,
A. P.
,
Sanger
,
J. M.
, and
Sanger
,
J. W.
,
2008
, “
Cardiac Myofibrillogenesis Inside Intact Embryonic Hearts
,”
Dev. Biol.
,
318
(
2
), pp.
236
246
.
35.
Torrent-Guasp
,
F.
,
Kocica
,
M. J.
,
Corno
,
A. F.
,
Komeda
,
M.
,
Carreras-Costa
,
F.
,
Flotats
,
A.
,
Cosin-Aguillar
,
J.
, and
Wen
,
H.
,
2005
, “
Towards New Understanding of the Heart Structure and Function
,”
Eur. J. Cardio-Thorac. Surg.
,
27
(
2
), pp.
191
201
.
36.
Knight
,
M. B.
,
Drew
,
N. K.
,
McCarthy
,
L. A.
, and
Grosberg
,
A.
,
2016
, “
Emergent Global Contractile Force in Cardiac Tissues
,”
Biophys. J.
,
110
(
7
), pp.
1615
1624
.
37.
Drew
,
N. K.
,
Eagleson
,
M. A.
,
Baldo
,
D. B.
, Jr.
,
Parker
,
K. K.
, and
Grosberg
,
A.
,
2015
, “
Metrics for Assessing Cytoskeletal Orientational Correlations and Consistency
,”
PLoS Comput. Biol.
,
11
(
4
), p.
e1004190
.
38.
Iyer
,
R. K.
,
Chiu
,
L. L. Y.
,
Reis
,
L. A.
, and
,
M.
,
2011
, “
Engineered Cardiac Tissues
,”
Curr. Opin. Biotechnol.
,
22
(
5
), pp.
706
714
.
39.
Chen
,
J. J.
,
Liu
,
W.
,
Zhang
,
H. Y.
,
Lacy
,
L.
,
Yang
,
X. X.
,
Song
,
S. K.
,
Wickline
,
S. A.
, and
Yu
,
X.
,
2005
, “
Regional Ventricular Wall Thickening Reflects Changes in Cardiac Fiber and Sheet Structure During Contraction: Quantification With Diffusion Tensor MRI
,”
Am. J. Physiol.
,
289
(
5
), pp.
H1898
H1907
.
40.
,
N. V.
,
Shashidhar
,
R.
, and
Chandrasekhar
,
S.
,
1971
, “
Orientational Order in Anisaldazine in the Nematic Phase
,”
Mol. Cryst. Liq. Cryst.
,
13
(
1
), pp.
61
67
.
41.
Pasqualini
,
F. S.
,
Sheehy
,
S. P.
,
Agarwal
,
A.
,
Aratyn-Schaus
,
Y.
, and
Parker
,
K. K.
,
2015
, “
Structural Phenotyping of Stem Cell-Derived Cardiomyocytes
,”
Stem Cell Rep.
,
4
(
3
), pp.
340
347
.
42.
Bone
,
P.
,
2008
, “
calcCircle
,” MathWorks, Natick, MA, accessed Mar. 2014–Mar. 2016, http://www.mathworks.com/matlabcentral/fileexchange/19083-calccircle
43.
Chernov
,
N.
,
2009
, “
Ellipse Fit (Direct method)
,” MathWorks, Natick, MA, accessed Mar. 2014–Mar. 2016, http://www.mathworks.com/matlabcentral/fileexchange/22684-ellipse-fit-direct-method/content/EllipseDirectFit.m
44.
Kroon
,
D.-J.
,
2010
, “
Snake: Active Contour
,” MathWorks, Natick, MA, accessed Mar. 2014–Mar. 2016, http://www.mathworks.com/matlabcentral/fileexchange/28149-snake—active-contour/content/InterpolateContourPoints2D.m
45.
D'Errico
,
J.
,
2012
, “
Inhull
,” MathWorks, Natick, MA, accessed Mar. 2014–Mar. 2016, http://www.mathworks.com/matlabcentral/fileexchange/10226-inhull
46.
Kuo
,
P.-L.
,
Lee
,
H.
,
Bray
,
M.-A.
,
Geisse
,
N. A.
,
Huang
,
Y.-T.
,
,
W. J.
,
Sheehy
,
S. P.
, and
Parker
,
K. K.
,
2012
, “
Myocyte Shape Regulates Lateral Registry of Sarcomeres and Contractility
,”
Am. J. Pathol.
,
181
(
6
), pp.
2030
2037
.
47.
Chiu
,
C.-L.
,
Digman
,
M. A.
, and
Gratton
,
E.
,
2013
, “
Cell Matrix Remodeling Ability Shown by Image Spatial Correlation
,”
J. Biophys.
,
2013
, p.
532030
.
48.
Feng
,
Z.
,
Qingming
,
H.
, and
Wen
,
G.
,
2006
, “
Image Matching by Normalized Cross-Correlation
,”
IEEE International Conference on Acoustics Speech and Signal Processing
(
ICASSP
), Toulouse, France, May 14–19, pp.
1977
1980
.
49.
Noorafshan
,
A.
,
Karbalay-Doust
,
S.
,
Khazraei
,
H.
,
Rafati
,
A.
, and
Mirkhani
,
H.
,
2014
, “
Spatial Arrangement of the Heart Structure: Application of Second-Order Stereology in Diabetic Rats
,”
Ann. Anat.
,
196
(
1
), pp.
20
25
.
50.
Hanson
,
B.
,
,
K.
,
Matsuura
,
K.
,
Robeson
,
S. M.
, and
Willmott
,
C. J.
,
1992
, “
Vector Correlation: Review, Exposition, and Geographic Application
,”
Ann. Assoc. Am. Geogr.
,
82
(
1
), pp.
103
116
.
51.
,
M. R.
,
Balachandran
,
K.
,
Capulli
,
A. K.
,
Golecki
,
H. M.
,
Agarwal
,
A.
,
Goss
,
J. A.
,
Kim
,
H.
,
Shin
,
K.
, and
Parker
,
K. K.
,
2014
, “
Engineering Hybrid Polymer-Protein Super-Aligned Nanofibers Via Rotary Jet Spinning
,”
Biomaterials
,
35
(
10
), pp.
3188
3197
.
52.
Deshpande
, V
. S.
,
Mrksich
,
M.
,
McMeeking
,
R. M.
, and
Evans
,
A. G.
,
2008
, “
A Bio-Mechanical Model for Coupling Cell Contractility With Focal Adhesion Formation
,”
J. Mech. Phys. Solids
,
56
(
4
), pp.
1484
1510
.
53.
Novak
, I
. L.
,
Slepchenko
,
B. M.
,
Mogilner
,
A.
, and
Loew
,
L. M.
,
2004
, “
Cooperativity Between Cell Contractility and Adhesion
,”
Phys. Rev. Lett.
,
93
(
26 Pt. 1
), p.
268109
.
54.
,
F.
,
Seidman
,
J. G.
, and
Seidman
,
C. E.
,
2005
, “
The Genetic Basis for Cardiac Remodeling
,”
Ann. Rev. Genomics Hum. Genet.
,
6
(
1
), pp.
185
216
.
55.
Graham
,
D. M.
, and
Burridge
,
K.
,
2016
, “
Mechanotransduction and Nuclear Function
,”
Curr. Opin. Cell Biol.
,
40
, pp.
98
105
.