Microscopic biophysical model of self-organization in tissue due to feedback between cell- and macroscopic-scale forces

We develop a microscopic biophysical model for self-organization and reshaping of artificial tissue, that is codriven by microscopic active forces between cells and an extracellular matrix (ECM), and macroscopic forces that develop within the tissue, finding close agreement with experiment. Microscopic active forces are stimulated by μm-scale interactions between cells and the ECM within which they exist, and when large numbers of cells act together these forces drive, and are affected by, macroscopic-scale self-organization and reshaping of tissues in a feedback loop. To understand this loop, there is a need to (1) construct microscopic biophysical models that can simulate these processes for the very large number of cells found in tissues, (2) validate and calibrate those models against experimental data, and (3) understand the active feedback between cells and the extracellular matrix, and its relationship to macroscopic self-organization and reshaping of tissue. Our microscopic biophysical model consists of a contractile network representing the ECM, that interacts with a large number of cells via dipole forces, to describe macroscopic self-organization and reshaping of tissue. We solve the model using simulated annealing, finding close agreement with experiments on artificial neural tissue. We discuss the calibration of model parameters. We conclude that feedback between microscopic cell-ECM dipole interactions and tissue-scale forces is a key factor in driving macroscopic self-organization and reshaping of tissue. We discuss the application of the biophysical model to the simulation and rational design of artificial tissues.

A key factor in the self-organization and reshaping of tissues is active feedback between microscopic forces generated by cells and macroscopic tensions in the extracellular matrix (ECM). Cells in tissues grow within ECM, a biopolymer network with viscoelastic properties. Cells generate active forces that locally distort ECM, and also respond to tensions in the ECM. We know from experiments [1], that even in the absence of initial tension, cells in tethered artificial tissues play an active role to remodel the ECM and generate tension. On the other hand, cells respond to this tension and align in a feedback loop. Thus, tissues are active materials that manifest self-organized order and spontaneous symmetry breaking due to the local reactions of cells to the properties of their environment and vice versa. If sufficient numbers of cells act in the same way, the resulting macroscopic tension reshapes the tissue, and also directs self-organization in the tissue, creating a microscopic-macroscopic feedback loop. Tension fields generated by these microscopic distortions and macroscopic reshaping of the ECM, lead to regions of both self-organization and disorder within engineered tissue samples [1][2][3][4].
The goal of this paper is to obtain greater insight into the processes of ECM-cell interaction, by developing a microscopic biophysical model that includes those interactions on a cellular scale, yet which is sufficiently simple to describe how these microscopic processes affect, and are affected by, self-organization on the macroscopic length scales of tissue. Therefore, we aim to develop a biophysical model that (1) is sufficiently complex to incorporate the physics of the subtle interplay between the extracellular matrix and active forces in cells, (2) can simulate the reshaping process and self-organization on a macroscopic scale, and (3) is sufficiently simple to describe large numbers of cells in tissue structures. The alignment of cells in tissues is similar to the nematic order found in liquid crystals, raising the possibility that microscopic models of the type used in statistical physics could be used to simulate artificial tissues. Our approach is different to previous work, with emphasis on a mixed phenomenological and biophysical model that is microscopic, yet fast to solve so that predictions can be made for the behavior of extremely large numbers of cells in tissues.
Continuum models have been developed to describe the alignment of stress fibers in artificial tissue constructs through both passive elasticity and the active forces generated by cells [5][6][7]. Legand et al. [8] have compared the continuum model of Refs. [5,6] with artificial tissue structures grown in microelectromechanical systems (MEMS), finding good agreement between model and experiment. Continuum models necessarily coarse grain cells and ECM into a single effective material containing the properties of both cells and ECM. Once this has been done, cells cannot be individually tracked, and cell alignments are assumed to follow stress fibers directly. A description of tissues as active continuum materials is well suited for modeling situations where the properties of individual cells are not important.
Microscopic models are complementary, but distinct to, continuum approaches. In a microscopic model, it is possible to model the complex self-organization of an ensemble of distinct cells in response to their environments, without the coarse graining associated with continuum models, and thus allowing simulations to be carried out that treat processes on cellular length scales. On the other hand, a typical challenge for microscopic models is carrying out simulations that are large enough to capture macroscopic scales.
To develop a microscopic model of ECM-cell interactions leading to self-organization, that can be simulated on macroscopic scales, we take inspiration from force-dipole models and contractile network models. A force-dipole approach could be used to simulate the interactions between cells and the ECM (see, e.g., Ref. [9]). Since it has been established that cells find an orientation consistent with the minimum work needed to obtain that configuration [10], it would then suffice to find the lowest-energy (ground) state of the theoretical model. However, the cells in such force-dipole models must be constrained to occupy nearly fixed positions in their substrates to avoid volume collapse relating to the form of the interaction (which has a 1/r 3 attractive form). Thus, such models cannot describe the reshaping of tissue that occurs simultaneously with self-organization. We also take inspiration from contractile network models, in which networks of springs can be used to model elastic effects in, e.g., biopolymers [11,12]. To date, the active forces due to cell-ECM interaction have not been introduced into contractile network models.
Our contractile network dipole orientation (CONDOR) model of cellular self-organization goes beyond previous work by combining force-dipole orientation models with contractile network models. Our microscopic model combines a complex network of elastic fibers that represents the ECM, with active forces that follow cell symmetries, to describe the active feedback between cells and ECM. This active feedback leads and responds to the macroscopic self-organization and reshaping of tissues. It is able to describe the tension fields, self-organization, and reshaping of tissues that arise from the active forces between cells and ECM. Spring networks are a popular way of representing elasticity in biopolymers [12], and force dipoles are a successful way of representing ordering of cells in tissues [9], but the two approaches have not yet been combined in this way.
Beyond the biophysical description of complexity in tissues, our CONDOR model could assist in the rational design of the scaffolds and molds that direct self-organization in engineered (or artificial) tissues. Artificial tissues are important in both pharmaceutical testing and regenerative medicine. In order to be useful for these applications, engineered tissues need to mimic the properties of real tissues. For example, it is critical that cells in artificial nervous and corneal tissues are self-organized so that they are highly aligned [13][14][15][16]. The development of new molds and scaffolds to encourage biologically relevant self-organization can take many months or years of painstaking trial and error. Thus, there is a need for new biophysical models capable of simulating the onset of self-organization in artificial tissues, which can be used for predicting which mold and scaffold designs have a good chance of success. Hence, our secondary goal is to generate techniques for detailed large-scale simulations of self-organization in artificial tissues that can be used as tools to assist the rational design of molds and scaffolds.
This paper is organized as follows. In Sec. II, the CON-DOR model and simulated annealing method for solving it are introduced. We present results from our microscopic biophysical model in Sec. III, which also includes comparisons between simulation results and existing experimental data. Finally, in Sec. IV we conclude and discuss the context of our microscopic biophysical model, the validation and calibration of our model, how our simulations could be used for the rational design of artificial tissue molds, and finally the outlook for extensions of the model.

II. METHOD AND CONDOR MODEL
In this section, we develop the CONDOR biophysical model that we use to describe how the active feedback between cells and ECM leads to reshaping of, and selforganization in, tissues, discuss how the parameters of the CONDOR model can be related to the cell density, and describe details of the simulated annealing method that we use to solve the CONDOR model.

A. Development of the CONDOR model for self-organization in engineered tissue
The biophysical model developed here is based on a threedimensional (3D) network of springs representing the ECM, with the interactions between cells and ECM introduced by changes in the undeformed length of springs dictated by the symmetry of the cells. The three-dimensional spring network has a face-centered-cubic (fcc) geometry. The cells in the artificial tissue structures that we simulate here are well separated, such that there is no direct interaction between cells; rather, the interaction is mediated via the springs representing ECM. The glial cells considered here (and many other cell types) are naturally elongated, and thus assumed to have approximately order 2 rotational symmetry. This is represented by an order 2 correction to the undeformed length of springs consistent with a force dipole. The network may, or may not, be subject to external stresses (depending on the system to be simulated). In all the cases considered here, the network is not subject to external stresses. A schematic of the model is shown in Fig. 1. A discussion of the form of the model begins by considering the energy of the contractile network representing the ECM in the absence of cells, which is given by a sum over Hooke's law for each individual filament, E = 1 2 i j κ i j (|l i j | − l u,i j ) 2 . Here, l i j is the displacement between cells at vertices i and j, l u,i j is the undeformed length of springs, and κ i j is the spring constant. In the absence of cells and external tension, the network is stress free. In all the examples considered in this paper, there is no external tension, so the network is stress free in the absence of cells. In the cases we consider, cells generate stress in the network. However, it is also possible to introduce stress by pulling on tethers that hold the tissue in place.
We introduce the effect of the active forces generated between cells and the network, via a modification of the undeformed length of the neighboring springs, l u,i j = l 0,i j − f ( 1,i j , 2,i j ). l 0,i j is the undeformed length of the bond in the absence of cells, and 1,i j and 2,i j are the angles between a bond and the cell orientations at the end of the bond. For simplicity, cells are positioned at the vertices between springs. Cells are able to remodel, pull on, and respond to the extracellular matrix [17], and as such the tissue sample, which is a material comprising both ECM and cells, is active and can generate its own stresses. It is noted that the modification we consider here represents the end result of the active processes, rather than the time-dependent details of such processes.
A form of f consistent with order 2 symmetry is where 0 < < 1 is a dimensionless constant controlling contraction perpendicular to the long axis of the cell. For each bond, cos( i j ) =l i j ·ŝ i , whereŝ i is a unitary vector representing the orientation of the cell andl = l i j /|l i j |. There are two length corrections, since there are cells at both ends of the bond (spring).
Thus, taking into account the effects of corrections to the length due to active forces generated by cells, the total energy of the CONDOR model is A schematic of the CONDOR model can be seen in Fig. 1. is a microscopic phenomenological parameter, which contains information about the level of contraction in the ECM caused by the activity of cells. It represents a drawing in of the ECM perpendicular to the orientation of the cell. Thus, it relates the active forces in the cells to the physical system of springs representing the ECM that connects the cells. We are not able to derive a form for this parameter, since it depends on cell-specific processes. In Sec. II B, we shall show how it is possible to establish a generic relation between the cell density and , by introducing a single free parameter that can be determined from comparison with experiment, such that predictions can be made.
The spring network used here is constructed on a regular lattice, and it is important to select a lattice with appropriate shear moduli to represent the ECM. Spring networks constructed from square or simple cubic lattices have zero shear modulus [12]. Thus our fully three-dimensional simulations are carried out using face-centered-cubic (fcc) lattices of cells, since the tetrahedral structures in these lattices ensure shear modulus. The shear and bulk moduli of the ECM, can be controlled more subtly through the ratio of next-nearest-neighbor (NNN) to near-neighbor (NN) coupling κ NNN /κ NN . Moreover, the inclusion of κ NNN improves the stability of the spring network against volume collapse.

B. Relating the CONDOR model to artificial tissue
In artificial tissue experiments, where cell-laden hydrogels are cast into molds, two key parameters control how cells interact with the ECM: (1) the cell density in the hydrogel and (2) the density of biopolymers in the hydrogel (for example, collagen). In this section, we describe how parameters of the CONDOR model and artificial tissue can be related.
Samples grown in cylindrical wells undergo larger contraction as seeding density increases, until contraction saturates at high densities [3]. To understand why contraction of artificial tissue samples saturates at large seeding densities, we consider how regions within which cells influence the ECM overlap as density increases. If density is low, the volume of influence 043217-3 increases linearly with cell density. However, at high density the volumes of influence around each cell start to overlap, so increasing the density further does not increase the overall volume of ECM that is influenced by the cells.
The influence of this overlap on the contraction can be estimated for a random distribution of cells, and used to arrive at a functional form to relate cell density ρ to . We denote the volume of influence for a single cell as V infl , and the volume of influence of all n cells as V tot,n . Note that V tot,n < nV infl , since the volumes of influence can overlap when cells are densely packed.
The mean additional volume of ECM, influenced by adding the nth cell, is V infl (1 − V tot,n−1 /V 0 ), where V 0 is the total volume of the tissue. This can be determined by considering the effect of adding an additional cell to a set of n − 1 cells. When the nth cell is added, the probability that its volume of influence overlaps with that of existing cells is V tot,n−1 /V 0 , which is the ratio of V tot,n−1 (the volume of influence of the existing n − 1 cells) and V 0 . Thus, the probability that the volume of influence of the new cell does not overlap is (1 − V tot,n−1 /V 0 ), and the mean additional volume of ECM influenced by the new cell can be calculated.
Thus, a recurrence relation is generated for the nth term in the expansion of the total volume of influence V tot,n , Repeated iteration then leads to a geometric series. Since R < 1, this series can be computed analytically for the volume of influence of N cells as So, assuming that ∝ V infl,tot , and since ρ = V 0 N = N/ρ 0 (we define ρ 0 = 1/V 0 , the density for a single cell within the tissue volume), Thus, the functional form relating the CONDOR model's contraction parameter to seeding density ρ is (noting that the maximum possible = 1) Without loss of generality, we define R = R ρ 0 /ρ 0 , where ρ 0 is a reference cell density. Since R and ρ 0 are independent, we fix R = 1/2 and allow ρ 0 , to set the density scale. So, With this choice, Eq. (8) has a linear form for ρ ρ 0 , saturates for ρ ρ 0 , and has the value = 1/2 when ρ = ρ 0 . Equation (8) may be inverted to determine ρ from :

C. Computational method
We determine cell orientations and positions using a simulated annealing approach to find the lowest-energy configuration of the cell-spring network. Bischofs and Schwarz have proposed that cells find an orientation consistent with the minimum work needed to obtain that configuration [10], i.e., the ground state. With a sufficiently slow anneal, the global minimum energy is guaranteed [18,19].
In the simulated annealing approach [18,19], a set of Monte Carlo updates are made, where in each update we modify either cell positions or orientations (or both), for either single cells or groups of cells. Cell positions and orientations are characterized using azimuthal and polar angles (θ and φ) and their coordinates (x, y, and z). Updates are accepted or rejected according to a Metropolis scheme, with probability exp(− E /η), where η is slowly reduced according to either an exponential or fast schedule. E is the change in energy due to an update.
During reorientation, both angles are updated either from a Lorentzian distribution [19] or a stepped distribution. For the Lorentzian distribution, angles are updated according to Here, r is a random variate with range [0,1), computed using the ran2 function from Press et al. [20]. The alternative step function update is Additional updates with θ = θ /10 and θ = θ /100 are also used. Similar updates are made for φ.
A major difference between this and previous work [10,21] is that cells can move to any location, thus allowing structural changes of the tissue to be explored. Lorentzian and step updates are made for x, y, and z. A fixed penalty is imposed on cells that move outside the working area of the mold. Within tethering regions, cells are not permitted to move, but can reorient.
Annealing is carried out using several anneal speeds, two anneal schedules and fixed/variable update ranges, with the lowest-energy configuration from the different attempts selected for analysis. Various annealing schedules were used, including fast annealing and exponential annealing. The initial η was selected by increasing η with fixed update range until a high acceptance rate was achieved. Various schemes for selecting the update range were used, including scaling of θ , etc., with η so that updates are accepted with roughly equal rate throughout the anneal. For a variable update range, the algorithm is considered to be finished when the ranges are less than a threshold, whereas for fixed range the algorithm finishes once the acceptance rate drops below a preselected value. Typical anneal times are a few CPU days.
Following the anneal, the lowest-energy configuration is used to calculate tensions in the springs, and contraction ratios. For consistency between contractile networks with different length scales, we introduce a dimensionless tension, T = T /κ NN l 0 . This measurement describes the tension in the environment of a cell.

D. Experimental data
We also carry out the task of validating our microscopic biophysical model, by comparing with experimental data that have been published elsewhere [3,4]. We have not performed  any new experiments for this paper, and emphasize that the primary purpose of this paper is to introduce a microscopic model of cell-ECM interactions in tissues. A summary of the raw experimental data from Refs. [3,4] that we have used in our comparisons can be found in the Supplemental Material [22].
Full experimental details have been reported previously, however, for convenience we provide a brief summary here [3,4]. Cells (C6 rat glioma cells or primary rat glial cells) were mixed with a neutralized solution containing acid-solubilized 2-mg/ml type I rat tail collagen and minimum essential medium then cast within multiwell plates (for contraction assays) or bespoke tethering molds (for alignment experiments). Cellular collagen gels were allowed to set for 10 min then immersed in culture medium and incubated for 24 h at 37 • C. During incubation, time natural cell-matrix interactions resulted in gel contraction and, in tethered gels, tensionmediated alignment of cells and matrix. Cellular gels were fixed and analyzed using immunofluorescence staining and confocal microscopy to visualize cell shape and position, with quantification of cellular alignment assessed using VOLOCITY software (Perkin Elmer, Waltham, MA).
There are two minor additional analyses carried out here. We convert the area contractions in Ref. [3] to a linear contraction for comparison with our model, by assuming that the resulting gels are circular and obtaining the radius of that circle, with error bars updated accordingly. For the tissue reported in Ref. [4], we reanalyze experimental cell orientations to provide information on cell orientations across the whole sample for comparison with theoretical predictions. This is done from the original confocal micrograph using VOLOCITY.

A. Untethered tissues in wells
In this section, we apply our microscopic model to one of the simplest 3D artificial tissue constructs, which has no tethers and can be said to be free floating. Free-floating tissue can be grown in circular wells (96-well plates). For convenient comparison with results from our microscopic model, Fig. 2 displays an image reproduced from Ref. [3], showing how the shape of the tissue is affected by the presence of cells. Initially, the whole well is filled with a mixture of collagen and cells. During the incubation period, the cells cause the tissue to contract. Cell-induced contractions in untethered artificial tissues increase with cell density.
Our computational simulations of self-organization in freefloating (i.e., untethered) artificial tissue in cylindrical wells, show increased contraction as gets bigger. Figure 3 shows top view theoretical calculations for the cell orientation in free-floating artificial glial tissue grown in 96-well plates for a variety of values. For the 96 well plates, the ratio of initial gel height to diameter is ∼0.1. Domains are visible in the cell orientations. Shapes are similar to the experiments shown in Fig. 2.
From these calculations, the contraction r/r 0 vs can be determined for free-floating gels in cylindrical wells, and is found to increase monotonically, saturating at ∼50% contraction. Here, r is the change in radius, and r 0 the initial radius. Figure 4 shows r/r 0 for simulations with various values for free-floating (i.e., untethered) artificial glial tissue grown in circular 24-well and 96-well plates [23]. The ratio of initial thickness to well diameter is ∼0.2 in the 24-well plates. Contractions in the horizontal plane are shown, since they correspond to measured experimental contractions from Ref. [3]. For free-floating gels, contraction tends towards ∼50% as → 1. The depth of the gel does not contribute strongly to the size of the contraction.
A linear interpolation of simulation results in Fig. 4, combined with experimentally measured contractions from Ref. [3], are used to determine how relates to ρ. Experimental data from Ref. [3] exist for two cell types, the C6 glioma line and primary rat astrocytes. Ranges of are determined from the experimental variance. Calculations corresponding to 96-well plates were used.
Using these values of ρ 0 to set the density scale, it can be seen that our microscopic model predicts contraction versus density curves that agree well with experiment (Fig. 6). Of course, this is to be expected, since a combination of theoretical and experimental values were used to find ρ 0 , however, agreement would not be possible if the functional form in Eq. (8) were not reasonable. Our theoretical predictions and the experimental measurements from Ref. [3] are both shown. We note that the area covered by the tissue is measured in Ref. [3], and we have converted this to an effective radius for comparison with our results. Fig. 4 combined with experimental data reported in Ref. [3] (labeled OR2015) to determine how relates to ρ, with ranges of determined from the contraction in 96-well plates, and the resulting plot is used to determine ρ 0 from using a least-squares fit of Eq. (8). The dashed line shows the fit and the dotted line the 95% confidence interval on the fit.

B. I-shaped tethers
In this section, we apply our biophysical model to the selforganization of engineered tissue in I-shaped clamps. Tethers can be used to guide tension in the engineered tissue. For   7. Image of tissues grown in an I-shaped mold, where reshaping of the ECM is modified by tethers (reproduced from Ref. [3]) for visual comparison with simulations from the theoretical model in Fig. 8. Cells are self-organized at the center of the sample (not shown). The contracted width of the tissue (highlighted with arrows) is reported in Ref. [3]. convenient comparison with our simulations, Fig. 7 shows an image of tissue grown in an I-shaped mold reproduced from Ref. [3]. The mold has clamping bars at the top and bottom of the image, within which there are multiple tethering points. In contrast to, for example, MEMS devices [8], the tethers that are simulated in this paper are passive, and do not exert any tension additional to that generated by cells. The biggest contraction of the tissue is towards the center of the sample, and there is also self-organization of cells in this region (not shown). The tether restricts contraction towards the top and bottom of the sample.
Tissue shapes predicted by our microscopic model agree well with experimental results from Ref. [3]. Figure 8 shows our theoretical predictions (shape and cell orientations) for artificial glial cells in an I-shaped clamp. These can be compared with experimental results (Fig. 7) from Ref. [3] (shape only). C6 cells of density ρ = 4 × 10 6 cells/ml form the artificial tissue.
Agreement is good between our theoretical predictions of contraction and experimental measurements from Ref. [3] (Fig. 9). Experimental data correspond to C6 cells tethered in an I-shaped clamp and infused with various concentrations of anti-CD49 antibodies. Anti-CD49 antibodies are used to block integrins, and therefore the ability for cells to bind to the ECM, leading to a significant decrease in the contraction of the tissue. If 96-well data (reported in Fig. 1 of Ref. [3]) are used in conjunction with Fig. 4 of this work to determine , then agreement between simulation and experiment with respect to contraction is excellent. We note that contraction of the hydrogel in the I-shaped clamp is smaller than predicted if Eq. (8) is used to set . The consequences of this for calibration of the CONDOR model will be discussed later. An important benefit of the microscopic model presented here is the ability to gain insight into quantities that would be difficult or impossible to determine experimentally, such as tension. In Fig. 10, we show theoretical predictions of the tension in free-floating artificial glial tissue grown in an I-shaped mold. In well-aligned regions, tension is low. Frustration associated with rapid changes of alignment with position leads to high tension close to the clamps. This high-tension region is Delta shaped. Insight into difficult-to-measure quantities could be used when designing bespoke clamps for artificial tissue with specific structures.

C. Pillar tethers
In this section, we discuss the application of our microscopic biophysical model to self-organization in molds where FIG. 9. Contractions of tissue in I-shape tethers from our theoretical model (reported here) show good agreement with experimental data (reported in Ref. [3]), if 96-well data are used to set . Points marked "experiment" correspond to the contraction of tissue tethered in an I-shaped clamp for hydrogels, made from C6 cells and infused with various concentrations of anti-CD49 antibodies. Agreement is poor if Eq. (8) is used to set . pillars tether the tissue. The relevant mold is a rectangular well with four pillars at each of the corners of the mold. As the hydrogel sets, the pillars hold the corners of the tissue in place, and direct tension as the tissue contracts under the influence of cells. Again, we note that the pillars are fixed and do not actively exert tension.
In a mold where pillars are used to tether the hydrogel, interactions between cells and ECM also lead to narrowing of the artificial tissue, as seen in confocal micrograph images [ Fig. 11(a), which is reproduced from Ref. [4]]. In Fig. 11(a), cells are stained red and the background with no cells is black. The artificial tissue starts as an untensioned rectangular type I collagen gel seeded with glial cells tethered by a set of four pillars, two to the left and two the the right of the image. These tethers show up as holes in the tissue sample after it is removed from the mold. It obtains a bowtie shape after 24 h of incubation. C6 cells of density ρ = 4 × 10 6 cells/ml form the tissue.
In order to determine the predictive power of our simulations, we determine cell orientations by analyzing the full confocal micrograph [ Fig. 11(b)]. The analysis highlights that the pillar-tethered hydrogel has a triangular (Delta) region just to the right of the two left-hand pillars, a region where cells are aligned along the short axis between the pillars and a region of glial cells aligned along the long axis towards the center of the sample. VOLOCITY (Perkin Elmer, Waltham, MA) is used to determine cell alignments from the image in Fig. 11(a). Yellow indicates cells aligned along the short axis, and blue indicates cell alignment along the long axis. Approximately 10 000 cells are identified in the confocal micrograph image. We note that the sample has been compressed under a cover slip.
Our computer simulations of cells in a pillarlike mold [ Fig. 11(c)] reproduce key features seen in the analyzed confocal image [ Fig. 11(b)]: Regions of self-organized cell alignment can be seen between the pillars and along the long axis at the center of the sample, and regions of mixed orientations are found where regions of different alignment meet. The positions of cells in the CONDOR model are fixed around the pillar, but permitted to change angle. All other cells are FIG. 11. Our simulations agree well with both the shape and cell orientation of artificial tissue tethered with vertical pillars reported in Ref. [4]. (a) A confocal micrograph of the engineered tissue reproduced from Ref. [4]. Tethering holes are visible on the left of the image and there is a characteristic bowtie shape. (b) Cell orientations extracted from the micrograph using VOLOCITY. (c) Theoretical results show self-organized cell alignment along the long axis and also along the short axis between tethers, with a region of various orientations at the interface between the two, as in the experimental data. (d) Moderate tension is found in the central region, which is aligned along the long axis. Regions of low tension are seen at the edge of the sample. Mixed tensions are found in the triangular region of mixed orientations. Regions of highest tension correspond to the region aligned along the short axis between the pillars. allowed to relax their positions. Results are shown for = 0.8 (which is consistent with the density of 4 × 10 6 cells/ml in the experimental data) and l 0 = 100 μm [ Fig. 11(c)]. Calculations contain 59 891 cells.

043217-8
Regions of highest dimensionless tensionT are associated with regions of rapidly changing cell orientation, as in the tissues simulated in I-shaped molds. The strongest tension can be found in the small region between the pillars, extending into a roughly triangular shape of higher tension along the long axis [ Fig. 11(d)]. So again, higher tension is associated with regions of rapidly changing alignment. The color represents the dimensionless valueT , with yellow representing regions of highest positive tension, red regions of moderate tension, and black regions of lowest tension. Contraction at the center of the sample is consistent between experiment and theory.

IV. DISCUSSION AND CONCLUSIONS
In this paper, we have introduced the CONDOR model, a microscopic biophysical model that combines contractile networks and force dipoles, to describe self-organization due to the feedback between cells and the extracellular matrix in tissues. In the CONDOR model, a network of springs represents ECM and the active forces generated by cells are represented by an additional contraction with a dipole symmetry. The lowest-energy configuration of the CONDOR model is found using simulated annealing. Results from our CONDOR model match experimental measurements of artificial tissue from Refs. [3,4] both on a microscopic level from cell alignments, and on a macroscopic level for the shapes of tethered tissues. We consider the ability to make macroscopic-scale predictions from a microscopic-scale model to be an important step forward for modeling artificial tissue.

A. Model validation and calibration
We have validated our microscopic biophysical model against existing experimental data that can be found in Refs. [3,4]. Our model performs well for a variety of situations, including free-floating gels in 96-well plates, gels clamped in I-shaped molds, and gels tethered using pillars. It is capable of good agreement with both the shape of the tissue and the orientations of cells within the tissue.
Calibration of the microscopic model requires the determination of ρ 0 , and is more challenging. One strategy is to fit the outputs of the model to data from free-floating gels, and we have done this using the data reported in Ref. [3]. Determination of from multiwell plates is particularly sensitive to statistical fluctuations in experimental results, due to the slow change in contraction with ρ, especially at larger ρ values. This indicates that very high accuracy experimental data are required to determine ρ 0 and thus calibrate the model. The experiments in Refs. [3,4] were carried out with the usual care to control confounding variables such as passage (which is an integer and can be matched exactly) and confluency, which is controlled to around 60%-70% [24]. However, even with this care, there can be large variations in tissue contraction, due not only to these variables, but also since the tissue can sometimes adhere to the sides of the mold in addition to the tethering points. This variation leads to sizable uncertainties in the value of ρ 0 for a particular cell type.
An alternative strategy for calibrating the microscopic model is to determine for specific batches using the contractions seen in multiwell plates. This strategy leads to very close agreement between theory and experiment if tethered molds are seeded with precisely the same cell batch.

B. Context
Our microscopic biophysical model combines contractile networks and force dipoles to predict self-organization and reshaping of tissue due to the microscopic active forces between cells and the ECM. We are unaware of other microscopic biophysical frameworks for modeling this process. However, we note similar approaches.
Pure force-dipole models do not incorporate the feedback between self-organization and reshaping driven by active forces that are found in our microscopic model (and real tissues). The symmetries of the active forces in our model are similar to force dipoles. Force dipoles have been a popular approach for describing the interactions between cells and the ECM [9,10,21,[25][26][27][28][29][30]. The force-dipole approach is sophisticated and has provided insight into the reasons that cells orient, especially in response to external stresses. Repulsive terms leading to finite compressibility of the viscoelastic ECM are not present. To stop collapse, cells in these theoretical models are constrained either to fixed positions, or small movements about such fixed positions, so reshaping of the tissue cannot be described.
With the exception of the microscopic model presented here, we are not aware of any contractile network models that have been modified to include the active forces between cells and ECM. In our microscopic model, the ECM is represented by a contractile network model. Contractile network models are common in biophysics, and can be used to represent many elastic media [12]. For example, contractile network models have be used to describe the shapes of both tethered single cells and bulk tissues [11,31,32].
In the context of our microscopic model, continuum mathematical models provide a complementary, yet distinct, approach for describing self-organization in tissues. We note continuum approaches to the alignment of stress fibers in artificial tissues, which can accurately predict the shapes of artificial tissue grown in MEMS systems [5,6,8]. These approaches are concerned with actin fibers and assume cell alignment along such fibers. We are not aware of any other continuum approaches that can describe the feedback between reshaping and self-organization in tissues.
Our CONDOR model has several key strengths, including the following: (1) an accurate description of the self-organization and reshaping of tissue that is driven by the feedback between cells and ECM; (2) a predictive capability for the shapes of tissues and the orientations of cells in those tissues; (3) simplicity allows for simulation of large numbers of cells, allowing macroscopic tissue sizes to be approached; (4) generalizability to different symmetry arrangements by adding higher-order multipolar corrections or other alternative symmetries to f (θ ) to allow the investigation of other biophysics; (5) extensibility to disorder, which can be achieved by including random bonds or by making corrections to the equilibrium distances of bonds; and (6) similarities to models common in statistical physics (such as the Heisenberg model of spin alignments, and models of liquid crystals), allowing the extensive machinery of statistical physics to be used to solve the microscopic model (e.g., Monte Carlo methods).

C. Potential applications of the model
Specific patterns of cellular self-organization are important in many tissues, making it challenging to grow artificial tissues that properly represent those in living organisms. The use of artificial or engineered tissue for regenerative medicine has the potential to significantly improve human health [2] and provide new biological models for drug testing [4]. The cells in real tissues are typically arranged in ways that may be difficult to reproduce artificially. For example, this is a particular problem for nerve tissue, which must be well ordered to function properly. In the nervous system, neurons are typically guided by glial cells, so well-ordered glial cells in artificial tissue could lead to well-directed neuron growth for, e.g., nerve repair. For such applications, favorable characteristics are high levels of cell alignment and cross section that is uniform along the length of the tissue, while regions of misalignment around tethering points should be minimized.
The design of artificial tissue with realistic characteristics is something that is very difficult and time consuming to do experimentally using trial and error, so the ability to predict outcomes would be extremely valuable. The CONDOR model described here is fast enough to assist the rational design of experimental molds and scaffolds that guide artificial tissue growth. The comparison with experiment carried out here indicates that for a design process, a broad range of should be used to obtain indicative results for specific mold designs.

D. Future prospects
Future work will increase sophistication of the CONDOR model to improve predictions of self-organization in tissue. There are several ways in which the model could be extended: (1) Improved cell response to tension in the ECM: The cells in the CONDOR model are always elongated (polarized), however, experimental systems have regions where cells are not elongated, particularly within the Delta region, so the contraction of the sample may be overestimated. The CONDOR model could be modified so cells are unpolarized in regions of low tension.
(2) Increased sophistication to improve the description of the ECM and analyze the role of disorder: For example, the polymers making up the ECM have many signatures of disorder. Randomized spring networks could therefore improve the description of collagen networks in the ECM. An intriguing possibility for further physics investigation is the existence of phase transitions between self-organization and disorder in the CONDOR model.
(3) Inclusion of time dependence: Cell representations in the model could be extended to include dynamic responses. This would allow the study of time-dependent ECM remodeling, the evolution of samples during their growth phase [3], or oscillatory strain stimuli (e.g., Refs. [33][34][35]).
Our microscopic model is amenable to all of these extensions and work is currently in progress. Overall, we expect that a wide range of interesting biophysics can be studied within extensions to the CONDOR model introduced here.