The structure/function relationship dogma is central to understanding how biological systems function. The idea is deceptively simple: understanding the structural organization of biological systems, from massive ecological systems to the shape of a single protein, reveals the function of the system. The concept is powerful enough to have inspired a 200+ year-long effort to describe the components of our biological universe in ever finer detail, beginning with the Linnean taxonomic system of cataloging organisms based on their structural similarities, and culminating with microscale descriptions such as the complete genomes of several organisms, including humans [1]. The reductionist approach to biological research has thus reigned supreme for generations, and as a result we now understand how the linear arrangement of nucleotides encodes the linear arrangement of amino acids, how proteins interact to form functional groups such as signal transduction and metabolic pathways, etc.

But at each level of biological organization, we reach a wall- having reduced the complex biological universe to a myriad of minute parts, we encounter new forms of complexity: data overload and the "curse of dimensionality [2]." Simply put, we've taken our biological machines apart but can't put them back together again- our ability to accumulate reductionist data has outstripped our ability to understand it. Thus, we encounter a gap in the structure/function relationship: having accumulated an extraordinary amount of detailed information about biological structures, we can't assemble it in a way that explains the correspondingly complex biological functions these structures perform.

This gap is especially evident at the level of tissues, where most diseases and injuries are manifest. Heart disease and cancer remain the top two causes of death in the United States. One fundamental characteristic of both diseases is tissue failure: namely, errors in the structural organization and function of cells in the affected tissues. Likewise, it is estimated that one in six US residents requires medical treatment for an injury each year [3], yet the process of wound healing is so complex it is difficult to accurately predict how quickly most serious wounds will heal [4,5]. Existing models of wound healing rely on clinically relevant, but somewhat superficial, measures of tissue state such as reduction in wound area, linear advancement of wound edge, pain, and ease of use [4-6]. In fact, despite a multitude of genetic screens, biochemical assays, and imaging techniques, the "gold standard" for diagnosis and evaluation remains the expert opinion of highly trained pathologists who scan samples of the tissues in histopathology slides. In other words, the human eye is currently the most accurate tool we have available for identifying telltale alterations in the structure and function of diseased and damaged tissues. And it is clear that human judgment is not fail-proof: thousands of diseases are misdiagnosed every year, costing hundreds of millions of dollars in wasted or ineffective medical treatment.

To improve diagnosis and treatment of diseases and wounds, we need a better understanding of how the tremendous numbers of cellular and subcellular parts are organized into functional tissues. One strategy for achieving this is to employ robust methods for describing complex systems, adapted from math and engineering disciplines far outside traditional biomedical fields [7]. Viewed from this perspective, tissue organization and function can be treated as a design optimization problem: what is the optimal arrangement of cellular constituents that achieves the best tissue performance?

When applied to problems of biological complexity, this design optimization approach is sometimes called systems biology. If one begins with the assertion that healthy, native tissues represent a design optimum, the task of systems biology is to identify the "control knobs" that govern tissue structure and function, and the specific "settings" of these knobs that yield an optimally functional (i.e., healthy) tissue. A second important assertion is that tissues are "self-correcting," in that when damaged, they are capable of generating an appropriate response that restores them to their optimal condition (i.e., wound healing): how are the control knobs "turned" to restore optimal function in a wounded tissue? All systems biology approaches therefore focus on defining three characteristics common to self-correcting systems: robustness (ability to maintain phenotypic stability in response to perturbation), modularity (clustering of components into functional "teams"), and, most importantly, emergent properties (behaviors unique to the entire system, and not found in any of its constituents)[8]. Systems biology approaches also share a concept known as iterative refinement, meaning that they cycle between perturbing a biological system, analyzing the data thus generated, and predicting how the system will respond to a new perturbation; these cycles learn relationships.

Since the late 1980s engineering design principles have been applied to living systems to create replacement tissues de novo [9]. In its most basic sense, an engineered tissue construct (ETC) is a three-dimensional assembly of one or more cell types suspended in an extracellular scaffold material and fed by soluble molecules, including growth factors, hormones, and nutrients [10]. Once assembled, the ETCs are intended to be implanted as replacements for damaged or diseased tissues. The results thus far have been promising [11], and some enjoy widespread use in the clinic [12].

Stem cells are very popular sources for the cellular component in these ETCs because they have the ability to proliferate (thereby populating the ETC with a high density of cells) and undergo differentiation once they reach their desired location and cell density. While embryonic stem cells retain the ability to differentiate into all tissues found in the adult, adult stem cells isolated from a number of different organs retain a more limited differentiation capacity. In either case, the promise is the same: stem cells offer the potential to define and manipulate fundamental principles of cell and tissue behavior, which in turn will uncover a new set of therapeutic targets for correcting errors in cell and tissue function [13].

Due to our very limited understanding of the principles governing tissue structure and function, most current tissue engineering typically follows a trial-and-error approach to design [14]. While persistence can pay off in the long run, the inefficiency of this approach remains one of the major barriers to widespread clinical application of ETCs [15]. Some stem cells used in ETCs also have the capacity to form tumors in vivo[16]. Until it is known how these cells decide which phenotype to adopt, this too presents a significant clinical challenge.

Systems biology approaches have been employed to help uncover the mechanisms governing differentiation and function of tissue stem cells [17-25] Comparatively little is known about the molecular control of human mesenchymal stem cells (hMSC). Originally described by Friedenstein [26], hMSC are a popular choice for musculoskeletal tissue engineering. hMSC are multipotent, self-renewing cells that can be isolated from the adult bone marrow [27]. They are capable of differentiating into at least three (osteogenic, chondrogenic, adipogenic) and perhaps as many as eight distinct lineages [28].

Due to its complex nature, unraveling stem cell differentiation requires a multi-stage approach. Reductionist studies have thus far identified a small number of potential regulators of this process [29-31], but fail to capture the global effects of these candidates on stem cell behavior. High-throughput, macro-scale studies (genomics, proteomics, etc.) are much better equipped to capture global changes in a complex system, and are the preferred choice for sampling dynamic changes in stem cell "state." But by themselves, these methods are typically not equipped to develop rigorous, testable hypotheses concerning the mechanisms governing this behavior [32].

The second stage of a systems biology approach to cell differentiation is to model the observed protein activity/gene expression changes as a function of the input stimuli. This is where so-called "traditional" biologists collaborate with experts in mathematics and computer science. Several modeling approaches have been applied. For example, Janes and Lauffenburger [33] illustrate how deterministic, probabilistic, and statistical learning models can be used to extract information about proteomic networks. The random nature of proteomic, genomic, and signaling data suggests that machine learning methods can also be used to model complex behavior such as protein-protein interactions [34,35]. For a compherensive review on this topic, we refer the reader to [20,36].

Faced with analyzing a wealth of data points as inputs, systems biologists turn to dimension reduction techniques using linear algebra. In linear algebra, a matrix is a two-way model used to describe linear relationships between the variables in two dimensions (e.g., rows and columns in a table). For example, in a gene expression experiment, the concentration of a chemical stimulant can serve as the row variables, and the resulting gene expression values can be the column variables. An entry in such a table would describe the gene expression level associated with a particular chemical stimulant concentration. However, it is now quite straightforward to generate data with three or more variables (also known as modes) (e.g., concentration of stimulant, protein phosphorylation level, gene expression level, duration of stimulant exposure, etc.) that cannot be represented by matrices.

Standard two-way dimension reduction techniques such as Singular Value Decomposition (SVD) [37,38] (see Figure 1), which organize data in a matrix form that incorporates time (e.g., concentration of stimulant × duration, protein phosphorylation × duration, or gene expression × duration), identify the linear relationships between modes in a pairwise fashion but cannot analyze three or more modes simultaneously. Uncovering meaningful patterns in a process as complex as cell differentiation, which obviously has far more than two modes, requires a more complex modeling approach.

This requirement can be satisfied by generalizing matrices to higher order models (e.g., moving from tables to n-dimensional cubes) to discover the multilinear relationships among data in datasets that have more than two different modes (i.e., multimodal data). Tensors are multidimensional arrays (also called n-dimensional cubes) ideally suited for multiway analysis of multimodal data. Figure 2 (top) illustrates the rows (X axis), columns (Y axis) and tubes (Z axis) of a sample of 3-way data array [39] while Figure 2 (bottom) illustrates the "slices" (known as elements) of the cube that correspond to the different modes in the data. In this way, one can generate several different matrices in the 3-way data cube. Though it is difficult to visualize, this approach can also be used to represent four, five, or more different modes simultaneously.

Having generated an "n-dimensional cube" data structure to represent the data, we now turn to the problem of how to fit a multiway model to the data and analyze the multilinear relationships between the modes. Two common models in multiway data analysis are Tucker3 [40-42] and PARAFAC [43].

Both methods model the original data by assembling a substantially smaller dataset representing the larger original data. While PARAFAC has not been applied to tackle the question of cell phenotype changes over time under different stimuli, Omberg et al. did use a Tucker3 (N-mode SVD) approach to model the time course of global gene expression in response to cell cycle inhibitors in yeast [44].

A higher-order tensor is a multiway dataset represented as T ∈ R^{I × J ×....N}, where M > 2. (Note that there are several notations used for representing multiway data sets such that both T and T∈RN1×N2×..NM refer to the same multiway array and will be used interchangeably.) For example consider a 3-way tensor i.e., a data cube that has locus links as its rows, gene ontology categories as its columns, and experimental conditions as its tubes. Below we will discuss three main techniques for analyzing such a tensor.

Matricizing (unfolding/flattening) rearranges three-way data as a matrix. Thus, it enables us to apply two-way analysis techniques, e.g., SVD, on a three-way dataset. As an example we illustrate matricization of a three-way tensor in the first mode (refer to Figure 3.). In the Tucker1 model the tensor is matricized in the mode of interest and SVD is performed on the corresponding matrix [45].

One of the most common multiway analysis techniques is the Tucker3 model. As depicted in Figure 4, a 3-way tensor T ∈ R^{I × J × K} is modeled as follows using a Tucker3 model :

T i j k = ∑ r = 1 R ∑ q = 1 Q ∑ p = 1 P G i j k A ip B jq C kr + E i j k

where P, Q and R indicate the number of components extracted from first, second and third mode (P ≤ I, Q ≤ J and R ≤ K), respectively. A ∈ R^{I × P}, B ∈ R^{J × Q} and C ∈ R^{K × R} are the component matrices. G ∈ R^{P × Q × R} is the core tensor and E ∈ R^{I × J × K} represents the error term. Different constraints such as nonnegativity, unimodality or orthogonality can be enforced on the component matrices. A Tucker3 model with orthogonality constraints on component matrices is a generalization of SVD from matrices to high-order datasets and is also called Higher-Order Singular Value Decomposition (HOSVD) [46] or multilinear SVD.

Tucker3 is the most flexible model among multiway analysis techniques. Although it suffers from core rotations, which result in non-unique solutions, a Tucker3 model may be preferred over other multiway analysis methods such as PARAFAC because of its flexibility in extracting a different number of components in each mode.

The simplest three-way model in terms of interpretation is Parallel Factor Analysis (PARAFAC) by Harshman [43] or Canonical Decomposition (CANDECOMP) by Carroll & Chang [47]. These two methods were originally proposed independently, but they employ the same model and are therefore considered completely equivalent.

An R-component PARAFAC model on a third-order tensor T ∈ R^{I × J × K} extracts three component matrices with R factors:

T i j k = ∑ r = 1 R A ir B jr C kr + E i j k

where R is the number of components extracted in each mode. A ∈ R^{I × R}, B ∈ R^{J × R}, and C ∈ R^{K × R} are the component matrices, and E ∈ R^{I × J × K} is the error term. An illustration of PARAFAC decomposition (Figure 5) gives more insight about how the model works. PARAFAC is a restricted Tucker3 model, where the core tensor is superdiagonal, meaning that only the entries G_iii in core tensor G can be nonzero for all i = 1..R.

The multiway modeling and analysis techniques presented above were applied to two systems biology problems: (i) discovering functional clusters of gene/protein expression during stem cell differentiation, and (ii) dynamics of hMSC osteogenic differentiation over time.

In [48] we identified 361 proteins (not expressed by the undifferentiated cells) that varied amongst the hMSC and hOST populations, and annotated the gene for each with its locus link number. We also mapped each gene in a gene ontology tree. The data are arranged as a tensor of type protein/gene locus link (LL) × category (Gene Ontology-GO) × stimulant, with dimensions (361 × 69 × 5), which we call the Proteomics Tensor (T). The stimulants were numbered as: collagen-I = sample 1, OS = sample 2, vitronectin = sample 3, hOST (positive control) = sample 4, hOST = sample 4, laminin-5 = sample 5.

We applied all three techniques reviewed above to model and analyze the T¹. Tucker1 model had (69, 26, 4) rank reductions for each mode, respectively to explain almost 90% variance of the data. These numbers were used to choose the core component numbers for our Tucker3 model. We also deployed a 12-component PARAFAC model based on the core consistency analysis [50].

The original proteomics data contained 71 gene ontology categories. We elected to drop the categories corresponding to Transmembrane Receptor Activity and Transmembrane Receptor Protein Tyrosine Kinase Signaling Protein because they contained very few genes. We have analyzed the data by using all three techniques considered in this work. Our analysis in Category mode also aims at identifying clusters of interest that may include categories that are nontrivially related.

We calculated a scatter plot of category names by projecting them on the first column vector of the category component matrix. The plotting is decluttered to provide a visualization of the clusters as shown in Figure 10. Similar to Tucker3 analysis, we plotted the category names by projecting them on the first column of the component matrix to obtain an intuition about the clustering by visualization as shown in Figure 11.

Expressing phenotypic changes over time requires multi-way data analysis, and thus is well suited to tensors. For example, Gaudet's "compendium" [55] includes signaling data triggered by different cytokines at various time points. One tensor in their compendium uses stimuli, measurements, and time as its modes. Here, we construct a 3-way tensor in which the rows are the gene IDs, the columns are the hMSC populations (negative control cells, hOST, and four samples each of the two stimuli), and the tubes are the replicates (all samples were repeated in triplicate). Over the course of the four time points we measured (days 1, 2, 4 and 5), we detected a total of 3153 genes with statistically significant (t-test, 0.05 confidence level) differential expression, relative to unstimulated hMSC. Thus the hMSC tensor T has 3153 × 10 × 3 elements as shown in Figure 15.

One objective of our analysis is to understand the evolution of the differentiation process over time, in particular the impact of different stimulants on this process. In Figure 16 we show the scatter plot of the data on the 1st component vector of the time & stimulant (population) mode. The 1st vector separates the two stimuli (static collagen I vs. stretched collagen I) and indicates that adding stretch accelerates osteogenic differentiation toward our target state, represented by fully differentiated hOST. This is consistent with our previous study which focused on a selected set of osteogenic marker genes and demonstrated that application of strain to collagen I induced a more rapid differentiation than static collagen I [28].

Our objective was to identify the outliers in the locus link mode and examine them in order to learn which genes are potentially important for the cell differentiation process.

For both Tucker3 and PARAFAC we computed the 98% concentration of the locus link numbers projected onto lower dimensions to capture the outliers in the remaining 2 percentile. Finally we took the intersection of the 2 percentile set of Tucker3 and PARAFAC.

When the data in this mode were projected onto the first component in our Tucker3 model (Figure 17), we observed a striking pattern: fully differentiated osteoblasts and unstimulated hMSC cultured on tissue culture plastic (TCP) lay at opposite ends of the graph, and the hMSC populations cultured on collagen in the presence or absence of mechanical strain were arranged in a temporal sequence, such that cells cultured for two days lie closer to the unstimulated hMSC, and cells cultured for five days lie much closer to osteoblasts.