
Vascular cells self-organize into unique structures guided by cell proliferation, migration, and/or differentiation from neighboring cells, mechanical factors, and/or soluble signals. However, the relative contribution of each of these factors remains unclear. Our objective was to develop a computational model to explore the different factors affecting the emerging micropatterns in 2D. This was accomplished by developing a stochastic on-lattice population-based model starting with vascular progenitor cells with the potential to proliferate, migrate, and/or differentiate into either endothelial cells or smooth muscle cells. The simulation results yielded patterns that were qualitatively and quantitatively consistent with experimental observations. Our results suggested that post-differentiation cell migration and proliferation when balanced could generate between 30–70% of each cell type enabling the formation of vascular patterns. Moreover, the cell-to-cell sensing could enhance the robustness of this patterning. These findings computationally supported that 2D patterning is mechanistically similar to current microfluidic platforms that take advantage of the migration-directed self-assembly of mature endothelial and mural cells to generate perfusable 3D vasculature in permissible hydrogel environments and suggest that stem or progenitor cells may not be fully necessary components in many tissue formations like those formed by vasculogenesis.
Citation: Jose E. Zamora Alvarado, Kara E. McCloskey, Ajay Gopinathan. Migration and proliferation drive the emergence of patterns in co-cultures of differentiating vascular progenitor cells[J]. Mathematical Biosciences and Engineering, 2024, 21(8): 6731-6757. doi: 10.3934/mbe.2024295
[1] | Floriane Lignet, Vincent Calvez, Emmanuel Grenier, Benjamin Ribba . A structural model of the VEGF signalling pathway: Emergence of robustness and redundancy properties. Mathematical Biosciences and Engineering, 2013, 10(1): 167-184. doi: 10.3934/mbe.2013.10.167 |
[2] | Peter Hinow, Philip Gerlee, Lisa J. McCawley, Vito Quaranta, Madalina Ciobanu, Shizhen Wang, Jason M. Graham, Bruce P. Ayati, Jonathan Claridge, Kristin R. Swanson, Mary Loveless, Alexander R. A. Anderson . A spatial model of tumor-host interaction: Application of chemotherapy. Mathematical Biosciences and Engineering, 2009, 6(3): 521-546. doi: 10.3934/mbe.2009.6.521 |
[3] | Marcello Delitala, Tommaso Lorenzi . Emergence of spatial patterns in a mathematical model for the co-culture dynamics of epithelial-like and mesenchymal-like cells. Mathematical Biosciences and Engineering, 2017, 14(1): 79-93. doi: 10.3934/mbe.2017006 |
[4] | Awatif Jahman Alqarni, Azmin Sham Rambely, Sana Abdulkream Alharbi, Ishak Hashim . Dynamic behavior and stabilization of brain cell reconstitution after stroke under the proliferation and differentiation processes for stem cells. Mathematical Biosciences and Engineering, 2021, 18(5): 6288-6304. doi: 10.3934/mbe.2021314 |
[5] | Avner Friedman, Yangjin Kim . Tumor cells proliferation and migration under the influence of their microenvironment. Mathematical Biosciences and Engineering, 2011, 8(2): 371-383. doi: 10.3934/mbe.2011.8.371 |
[6] | Wing-Cheong Lo, Ching-Shan Chou, Kimberly K. Gokoffski, Frederic Y.-M. Wan, Arthur D. Lander, Anne L. Calof, Qing Nie . Feedback regulation in multistage cell lineages. Mathematical Biosciences and Engineering, 2009, 6(1): 59-82. doi: 10.3934/mbe.2009.6.59 |
[7] | J. Ignacio Tello . On a mathematical model of tumor growth based on cancer stem cells. Mathematical Biosciences and Engineering, 2013, 10(1): 263-278. doi: 10.3934/mbe.2013.10.263 |
[8] | Elena Izquierdo-Kulich, Margarita Amigó de Quesada, Carlos Manuel Pérez-Amor, Magda Lopes Texeira, José Manuel Nieto-Villar . The dynamics of tumor growth and cells pattern morphology. Mathematical Biosciences and Engineering, 2009, 6(3): 547-559. doi: 10.3934/mbe.2009.6.547 |
[9] | Fiona R. Macfarlane, Mark A. J. Chaplain, Tommaso Lorenzi . A hybrid discrete-continuum approach to model Turing pattern formation. Mathematical Biosciences and Engineering, 2020, 17(6): 7442-7479. doi: 10.3934/mbe.2020381 |
[10] | Samantha L Elliott, Emek Kose, Allison L Lewis, Anna E Steinfeld, Elizabeth A Zollinger . Modeling the stem cell hypothesis: Investigating the effects of cancer stem cells and TGF−β on tumor growth. Mathematical Biosciences and Engineering, 2019, 16(6): 7177-7194. doi: 10.3934/mbe.2019360 |
Vascular cells self-organize into unique structures guided by cell proliferation, migration, and/or differentiation from neighboring cells, mechanical factors, and/or soluble signals. However, the relative contribution of each of these factors remains unclear. Our objective was to develop a computational model to explore the different factors affecting the emerging micropatterns in 2D. This was accomplished by developing a stochastic on-lattice population-based model starting with vascular progenitor cells with the potential to proliferate, migrate, and/or differentiate into either endothelial cells or smooth muscle cells. The simulation results yielded patterns that were qualitatively and quantitatively consistent with experimental observations. Our results suggested that post-differentiation cell migration and proliferation when balanced could generate between 30–70% of each cell type enabling the formation of vascular patterns. Moreover, the cell-to-cell sensing could enhance the robustness of this patterning. These findings computationally supported that 2D patterning is mechanistically similar to current microfluidic platforms that take advantage of the migration-directed self-assembly of mature endothelial and mural cells to generate perfusable 3D vasculature in permissible hydrogel environments and suggest that stem or progenitor cells may not be fully necessary components in many tissue formations like those formed by vasculogenesis.
A major obstacle in the development of tissue engineered products for clinical applications is the challenge of generating perfusable vasculature, in vitro. This issue is amplified in both importance and scale when building larger or highly complex organs. While the decellularization of adult organs such as heart [1,2], lung [3,4,5], kidney [6,7], and liver [8,9] retains the highly branched vascular architecture, methods for cell seeding preformed matrix at physiologically-required cellular densities remains a challenge. Another promising approach to generating perfusable vasculature is seeds mature vascular cells into 3D hydrogels. Several groups have developed microfluidic platforms that permit the self-assembly of endothelial cells into perfusable vasculature [10,11,12] and direct anastomosis of that vasculature within the platforms microstructures. In these studies, co-cultures of human umbilical vein endothelial cells (HUVECs) and normal human lung fibroblasts (NHLFs) or pericytes (PCs) are seeded as single cells into fibrin or collagen gels [13,14,15]. Over 4–7 days, the migration and self-assembly of these cells leads to the formation of perfusable vessels. A third approach to generating perfusable vasculature uses developmental or differentiation methods enabling endothelial cell (EC) [16,17] and smooth muscle cell (SMC) [18,19,20] fate within organoid-like structures supplemented with growth factors, specialized cell culture mediums, and/or mechanical signaling [21] to guide the formation and development of perfusable vessels [22].
Our laboratory's chemically-defined differentiation protocols have proven highly effective in deriving mouse embryonic stem cells (ESCs) and human ESCs and induced pluripotent stem cells (iPSCs) [23,24,25,26] into vascular progenitor cells (VPC), ECs, and SMCs (Figure 1A). While obtaining purified differentiated vascular cells is an important first step, the ultimate goal is to simultaneously direct the emergence of vascular structures from differentiating cells as seen during normal development. Indeed, VPC outgrowths have been observed to self-organize into micropatterns (Figure 1B) with EC clusters loosely surrounded by SMCs [27]. However, in order to direct vascular patterning in 3D, one must understand and control the interplay of various processes involved in cell organization: Proliferation, migration, and differentiation. While the details of such strategies are highly system-dependent given the complexity and inter-connectedness of the processes, we sought to first understand the relative importance of these processes in patterning. For example, do cell patterns emerge from directed differentiation from neighbor cells or the migration and rearrangement of post-differentiated cells (Figure 1C)? Here, we examine the relative importance of physiologically relevant range of parameters using a computational model motivated and calibrated by our experimental cultures of mouse ESCs in 2D [27].
Several computational models exist for simulating the self-organization of multicellular tissues [28], including on-lattice and off-lattice models that account for cell adhesion, proliferation, and short or long range signaling. Depending on the specific question, different models exhibit distinct advantages. For example, vertex models are excellent for exploring adhesion, proliferation, cellular forces, and cellular geometries, but become computationally expensive with more cells [28]. Likewise, finite element models (FIE) have been used in modeling environments where cell geometries and cellular forces are important [29,30]. Many other models have explored cell differentiation [31,32], migration [33], proliferation [34], and combinations of these parameters [35,36,37,38,39,40,41]. Here, we chose an on-lattice, stochastic, population-based model that uses ordinary differential equations (ODE) to represent the spatial-temporal dynamics of cell densities (# of cells per lattice site) evolving through cell proliferation, migration, death, differentiation, and cell-to-cell signaling. Such an approach allows effective parametrization of the important processes, while also monitoring the spatial structure of the cellular populations over time. A coarse-grained approach also avoids replicating dynamics at the single cell level increasing computational efficiency in simulating the large cell numbers (104–105) observed in our experimental cell cultures.
Here, we simulated the time evolution of the system for a wide range of physiologically plausible parameters and output reporting the fractions and spatial distributions of the emerging vascular cells (ECs or SMCs). By analyzing simulations within physiologically relevant domains calibrated with experimental data, we extract practical information from each of the simulated parameter sweeps. The results suggest that the vascular-like spatial patterns (defined as EC clusters surrounded by SMCs) emerge when the fraction of differentiated ECs lies within a "zone of co-emergence" with well-balanced numbers of both cell types. Moreover, for physiologically relevant ranges of proliferation and migration, the distribution of simulated EC cluster diameters was consistent with experimental observations. For the parameter ranges relevant to our system, we found that the spatial distributions of different cells are more sensitive to differences in proliferation and migration rates between cell types compared with differences in intrinsic or induced differentiation rates. Consequently, it is the proliferation and migration rates that appear to be most responsible for the establishment of the observed micropatterns within our differentiating co-cultures. These results support the current self-assembly vasculogenesis observed in mature post-differentiation cell co-cultures when seeded withing hydrogels [10,11,12] and may aid the rational design of co-developing vasculature within 3D tissues.
The mouse embryonic stem cell (mESC) lines used for these studies included mESC-R1 (ATCC). The mESC-R1 were cultured on 0.5% gelatin in serum-free medium containing Knockout Dulbecco's Modified Eagle Medium (KO-DMEM; Invitrogen), 15% Knockout Serum Replacer (KSR; Invitrogen), 50 units/mL and 50 ug/mL of Penicillin-Streptomycin (Invitrogen), 1Χ Nonessential Amino Acids (Invitrogen), 2 mM L-glutamine (Invitrogen), 0.1 mM 2-mercaptoethanol (Calbiochem), 2000 Units/ml of leukemia inhibitory factor (LIF-ESGRO; Chemicon), and 10 ng/ml of bone morphogenetic protein-4 (BMP-4; R & D Systems). Full media changes occurred every other day and cells were passaged every four to five days and reseeded at a density of 104 cell/cm2.
Mouse ESC-R1 (ATCC) were differentiated into EC and SMC using our laboratories two staged serum-free induction protocols [25]. Briefly, mESC were induced on 50 μg/mL fibronectin (Corning) coated plates (BD Biosciences) under our stage 1 induction medium containing alpha-MEM (Cellgro), 20% knockout serum replacement (ThermoFisher), 50 units/mL and 50 ug/mL of penicillin-streptomycin (ThermoFisher), 1Χ nonessential amino acids (ThermoFisher), 2 mM L-glutamine (ThermoFisher), 0.05 mM 2-mercaptoethanol (Calbiochem), 5 ng/mL BMP-4 (Peprotech), and 30 ng/mL of VEGF (Peprotech). After 2 days in culture, vascular progenitor cells (VPCs) were live stained for Flk-1 expression (1: 200, Biolegend) for 30 minutes, washed, and sorted by Fluorescence-activated cell sorting (BD FACS Aria Ⅲ). Positive cells were replated on 50 μg/mL fibronectin for an additional 4 days under stage 2 specific differentiation medium consisting of 70% alpha-MEM (Mediatech), 30% DMEM (Invitrogen), 2Χ Nutridoma CS (Roche), 50 units/mL and 50 ug/mL of penicillin-streptomycin (Invitrogen), 1Χ nonessential amino acids (Invitrogen), 2 mM L-glutamine (Invitrogen), 0.05 mM 2-mercaptoethanol (Calbiochem), and supplemented with 50 ng/mL of basic fibroblast growth factor (bFGF) as previously optimized for mESC-R1 induction into ECs [25].
Visualization of the cell micropatterning was conducted by immunofluorescent staining of VPCs outgrowths on day 4 post purification (Figure 1B). Briefly, cells were fixed with 4% paraformaldehyde (Tousimis) and permeabilized with 0.5% Triton X-100 (MP Biomedicals). Nonspecific binding was prevented using 1% bovine serum albumin (Sigma). Conjugated CD31 PE, EC stain (BD Biosciences) and primary antibody CNN1 (SMC stain, Sigma) were added and allowed to stain overnight at 4˚C. Cells were rinsed before addition of secondary antibody, Alexa Flour 488 (Thermofisher), and DAPI. Cells were incubated for an additional hour before final rinse and imaging via fluorescence microscopy (Nikon TE2000-U). For F-actin staining we used an Alex Fluor 488 phalloidin stain (Invitrogen).
The diameters of EC clusters were calculated by a custom MATLAB script. Briefly, imported images, either experimentally obtained images of DAPI stained cells or computationally simulated clusters, were smoothed via medfilt2 function, and then turned into binary (imbinarize) images via Otsu (graythresh) thresholding. An outline of the clusters was then reconstructed via subtraction of the dilated and eroded binary images. The regionprops function was then used to record the Filled Area of the clusters (Figure S1). Comparison between experimental and simulated clusters, was done by calculating an effective cluster diameter, D = 2*(Filled Area /π)0.5, for both conditions.
Our computational model considers the spatial and temporal distributions of three different cell types within a simulated square lattice. We denote by XA, XB, and XC the normalized cell density (e.g. XA = NA/Nmax, where NA is the number of A cells and Nmax is the minimum number of cells of all types put together allowed at each lattice site leading to X belonging to the interval [0, 1]) of VPCs, ECs and SMCs, respectively. Their dynamics are governed by the cell's migration, their rates of proliferation, differentiation, and cell death (Figure 2A–C). In our model, each cell can migrate between neighboring lattice sites at a preset migration value, Jθ (here θ represents either A, B or C corresponding to cell identities of VPC, EC, and SMC, accordingly). Each cell type can also proliferate (double), and die at a specific rate, δθ and μθ, respectively. Additionally, VPCs can differentiate into either ECs or SMCs at rates of αB and αC, respectively. The cellular dynamics are then a combination of a cell's migration, these rates, and stochastic noise, Sθ, leading to a set of three ordinary differential equations (ODEs), Eqs (1)–(3), one for each of the three changing cell populations (XA, XB, and XC).
˙XA=((δA+S)(1−∑θ∈A,B,C[Xiθ])−μA)XA−(αB+αC)XA+FA | (1) |
˙XB=((δB+S)(1−∑θ∈A,B,C[Xiθ])−μB)XB+αBXA+FB | (2) |
˙XC=((δC+S)(1−∑θ∈A,B,C[Xiθ])−μC)XC+αCXA+FC | (3) |
here, noise is captured by a stochastic addition to the combined rates of proliferation and death. The stochastic noise, Sθ, is drawn from a normal distribution and is assumed, for simplicity, to be the same for all cell types. The factor of (1−∑θ∈A,B,C[Xiθ]) modulates the growth rate relative to available space and accounts for the carrying capacity of each lattice site i.e., the sum of the normalized cell densities cannot exceed 1 and the growth rate vanishes as the total density approaches 1. The last term in each equation accounts for the flux of cells, Fθ, into a lattice site which depends on the cell's migration rate, Jθ. In the absence of any biases or interactions, these fluxes are proportional to the cell density gradient for a specific cell type, summed over all nearest neighbor sites, with a constant of proportionality reflecting the migration of the cell type Eq (4).
Fiθ=∑j∈n.nofiJθ(Xjθ−Xiθ) | (4) |
The range of physiologically relevant migration, proliferation, and differentiation rates we used were derived from various literature sources. For example, prior studies have reported that ECs and SMCs can have a maximum displacement distance of 18 and 44 μm over the course of an hour, respectively [33,42]. Additionally, proliferation rates for mouse ECs have been reported to have a minimum doubling time of 19 hrs [43], while 22 hrs has been reported for rat aortic vascular SMCs [44]. Last, the rate of differentiation can be considered as the length of time for specific EC or SMC marker expression. Some studies have reported initial EC marker expression at as early as 24 hrs post induction [45] and 36 hrs for SMC marker expression [46], while most studies do not see differentiation makers emerging before 5–7 days [23] from mouse stem cells and 10 to 20 days from human stem cells [47]. Using these physiologically relevant values for our parameters (see Table 1; Supplementary video V1), our simulations were evolved in time using the Euler (first order Runge-Kutta) method implemented in Python [48]. Each cell's population density, at each lattice site, is tracked over time (recorded every 1hr in time) and can be reconstructed to create a 2D spatial representation of the evolving cell populations using MATLAB [49] (Supplementary video V2 and Figure S2A–C).
Variable | Denotation | Experimentally obtained values | Corresponding simulation value(s) | Parameter sweep | References |
Diffusion constant | Jθ | 14 μm2/hr | 0.0079 | 0–0.1; 0.001 step size | Supplemental video V1 (37–38) |
Proliferation rate | δθ | ~40 hr | 0.025 | 0–0.1; 0.001 step size | (39–41) |
Differentiation rate | αθ | ~62.5 hr | 0.016 | 0–0.1; 0.001 step size | (19, 20) |
Stochastic noise | Sθ | A loss or gain of up to 1 cell per lattice site is incorporated every 10 hrs | ± 0.01 | N/A | NA |
Paracrine signal strength | β | amplification of differentiation rates by 0.5–1.5 | ± 0.5, ± 0.016, 0 | N/A | NA |
Under the simplest assumption, the various rates (in Table 1) for specific cell types are fixed. However, these rates can also be modulated by different types of external cell interactions, such as contact inhibition and paracrine signaling [25,50]. To explore how these external processes affect patterning, we accounted for their promotion or suppression of the corresponding rate constants for migration, proliferation, and differentiation.
To model the effects of contact inhibition, we assumed that a cell's migration rate would simply decrease linearly with the number of cells in its immediate neighborhood (summed over the site and all nearest neighbor sites) [51]. Within this scheme, we examined two distinct possibilities, specific and nonspecific cell adhesions, as possible regulators (Figure 2D). For specific homotypic cell adhesions, we assumed that cells can only sense and interact with other cells of the same cell identity, leading to a modified migration rate,
Jθ=Joθ(1−(Xiθ+∑j∈n.nofiXjθz)) | (5) |
where i, j are site indices and z is the coordination number of nearest neighbor lattice sites plus the current site (z = 9). The migration rate thus decreases from a maximum for isolated cells to zero when the cell density of the same type is maximal. For nonspecific cell adhesions, we assume cells can sense all other cells in their local neighborhood (leading to Eq (6)).
Jθ=Joθ(1−(∑θ∈A,B,C[Xiθ+∑j∈n.nofiXjθ]z)) | (6) |
here, the migration rate decreases from a maximum for isolated cells to zero when the overall cell density for all cells is maximal. Combinations with different cell types employing different sensing mechanisms were also explored. For example, ECs could be modeled to nonspecifically sense their neighboring cells, leading to greater contact inhibition, while SMCs could be assigned to specifically sense other SMCs (consistent with homotypic sensing). It is to be noted that patterning arising from variations in adhesion is consistent with the differential adhesion hypothesis [52], which states that cells with strong attractions will cluster closely together while those with weaker attractions will surround them. Here we account for the migration of cells being differentially affected by differences in cell adhesion dynamics mediated by specific vs nonspecific sensing. An example of a strong (homotypic) cell-cell adhesion would be vascular endothelial cadherins (VE-cadherins) that hold ECs together [53] with a binding force between 35–55 pN [54], while E-cadherin which can binds ECs to SMCs is an example of a nonspecific cell adhesion [55] with a binding force between 32–48 pN [56].
A similar approach was implemented for modifying the proliferation rates based on homotypic and nonspecific cell adhesion. We assumed that proliferation would linearly decrease as the number of neighboring cells increased, in accordance with evidence demonstrating that the mitotic rate may be arrested when a critical cell confluency is achieved. This type of contact inhibition emerges from an increase in mechanical interactions resulting from the loss of available space, thus constraining cellular dynamics including cell division [41]. It is also possible for arrest to occur in a specific manner such as being mediated by the EC specific Notch signaling pathway that inhibits tip-EC proliferation during angiogenesis [57]. Therefore, we modified the proliferation rates appropriately to account for two different sensing mechanisms: 1) Specific, homotypic, cell sensing:
δθ=δoθ(1−(Xiθ+∑j∈n.nofiXjθz)) | (7) |
and 2) nonspecific cell sensing:
δθ=δoθ(1−(∑θ∈A,B,C[Xiθ+∑j∈n.nofiXjθ]z)) | (8) |
These proliferation rate dependencies are analogous to Eqs (5) and (6) for the migration rates, with a coordination number (z) equal to 9.
For differentiation, we examined how chemical signals and growth factors produced by the cells influence VPC fate decisions, modeling vascular paracrine signaling [58]. Such signaling is observed in EC and SMC recruitment during early vascular development, with ECs secreting platelet-derived growth factor-b when recruiting SMCs [53]. We explored two possible mechanisms for guiding VPC differentiation into ECs and SMCs. The first mechanism, designated same cell-directed differentiation, involves a committed cell, EC or SMC, inducing a neighboring uncommitted VPC to differentiate into a cell of the same identity. The second mechanism, designated alternate cell-directed differentiation, assumes a committed cell, such as an EC, influencing a neighboring VPC to instead differentiate into the alternative cell type, in this example inducing it into a SMC (Figure 2E). For simplicity we assumed an isotropic distribution of paracrine signals for both cell types, where the diffusion range is equivalent to one lattice unit. Given that our lattice size is 79 μm (see next section), this is a reasonable upper bound for diffusive signal propagation. The modified differentiation rate is then given by
αθ=αoθ(1±β(Xiθ'+∑j∈n.nofiXjθ'z)) | (9) |
where the choice of EC or SMC, for θ and θ', determines the type of cell-directed differentiation mechanism (Figure S3) and β represents the strength of the signaling from cells in neighboring lattice sites.
At time 0 (t = 0), a fraction of our 50x50 lattice is randomly seeded with VPC population densities, XA, with values taken from a normal distribution between [0, 1]. Here each unique lattice site supports up to 10 cells and has a lattice unit size of 79 μm (Figure S2A). This choice of lattice size was motivated by the average cell size of our VPCs, measured to be 25 ± 7 μm across (Supplementary video S1, Figure S2D–E), whereby each lattice site is significantly larger than 1 cell diameter, so that our population-based approach remains valid, and yet small enough to visualize the population driven micropattern features consisting of mature cell clusters with effective diameters of 340 ± 110 μm (Figure S1G). The fraction of sites seeded was chosen such that the seeded density for the array was consistent with experimental cell seeding densities of VPCs (of 104 cells/cm2) [25]. A migration constant value was derived from the average distance traveled by randomly migrating human umbilical vein ECs, of about 14 μm over an hour [33]. Taking it to be a random walk process, we can compute an effective diffusion constant, D, from mean-squared displacement (<r2 > ~(14 μm)2) and the time taken (t), as < r2 > /4t = D. Non-dimensionalizing D using the lattice site dimensions (79 μm) as the unit of length and 1hr as the unit of time, we obtained a corresponding dimensionless migration constant of Jθ = 0.0079. Additionally, we chose to fix the baseline (unmodified) values for proliferation and differentiation rates within physiological ranges (Table 1). These rates were calculated from the typical times associated with the given processes. We set the rate of proliferation for all cells to 40 hrs [43,44,59], corresponding to a dimensionless simulation value of δθ = 0.025. The differentiation rate was set to 62.5 hrs [27] thus falling between days 2–3 of post stage 2 differentiation, corresponding to a simulation value of αθ = 0.016. Once the simulation is initialized, the three first order ODEs, Eqs (1)–(3), are evolved in time using the Runge-Kutta method, and Xθ values are recorded at every lattice site at 1-hr intervals for the duration of the simulation (Figure S2B, C). At t = 96 hrs, corresponding to experimental data at day 4, the recorded information is quantified and analyzed. We recorded the populations size for each cell type, their spatial distribution, and the fraction of ECs and SMCs at t = 96 hrs, under each specific condition. Supplementary video V2, shows a simulation evolving in time, for a particular set of parameter values. We conducted parameter sweeps over regions of parameter space that include experimentally relevant and physiologically plausible (see Table 1) values for each set of parameters. Simulations were run a total of 10 times for each set of parameters for statistical analysis.
Using our stochastic on-lattice population-based model, we will explore the role that migration, proliferation, and differentiation rates have on the type of patterns that emerge from VPC (A cells) as they differentiate into ECs (B cells) and SMCs (C cells). We will also explore the cluster dynamics of a single EC cluster under various migration rates, and end by simulating different types of cell-to-cell signaling such as contact-inhibition and paracrine signaling.
We first explored the role that cell migration, for EC (B cells) and SMC (C cells), has on the emerging micropattern. Our parameter sweep explored physiologically relevant migration values, between 0–18 μm2/hr for ECs and 0–44 μm2/hr for SMCs [33,42] corresponding to in-silico values of JB = 0–0.013 and JC = 0–0.078, respectively. We varied the ECs' (JB) and SMCs' (JC) migration rates while holding the VPC's (A cells) migration (0.0079) constant and setting all other parameters, such as differentiation (0.016), proliferation (0.025), and noise (0.01) to be the same for all three cell types (see Table 1). The fractions of ECs at t = 96 hrs were then plotted over this parameter space (Figure 3A). We verified that roughly equal migration values for JB and JC, generated relatively equal population densities of ECs and SMCs (EC fraction~0.5). The spatial distribution of cells in this regime also revealed micropatterning, defined as spatially separated regions dominated by one cell type or the other (Figure 3B i-iv). We observed that the micropatterning persists but becomes less distinct as the overall migration rates increased (compare Figure 3B i to iv). Moreover, micropatterning was observed over a broad range of EC fractions, roughly between 0.3 to 0.7, demarcating a "zone of co-emergence" (Figure 3A).
This zone of co-emergence can be quantified by cell type asymmetry (Figure S4A), where asymmetry is defined as (BT - CT)2/(BT + CT)2 where BT and CT are the total number of ECs and SMCs in the simulated space at t = 96 hrs. Here values of asymmetry close to 0 indicate similar population densities of ECs to SMC, while values closer to 1 denote the dominance of one cell type. The corresponding zone of co-emergence can be distinctly observed at asymmetry values up to 0.3 (Figure S4A). Interestingly, in our laboratory's experimental studies on these same differentiating co-cultures of VPCs [27], we reported EC fractions to be 0.30 at day 3 and 0.25 at day 7 with distinct patterns being observed at day 7. The standard deviations from the stochastic effects introduced are averaged over 10 simulations and were observed to be insensitive to repeated runs, with any variance localized within the zone of co-emergence (Figure S4B). Next, we looked at the cases where one cell type migrated an order of magnitude faster than the other cell type. Not surprisingly, we saw that the faster cell type dominates and outcompetes the slower cell type for space (Figure 3B, v and vi). However, even a more modest difference of about 50% in the migration rates allowed the faster cell type to dominate (Figure 3B vii and viii). In fact, upon further investigation, the zone of co-emergence, is highly sensitive to differences in migration rates and with differences as small as roughly half a cell size per hour (~4.5 μm/hr), we observe a significant dominance of the faster cells (Figure S5).
The region of parameter space within the zone of co-emergence where micropatterning occurs was next explored. To do this, we quantified the spatial separation between the cell types, by measuring EC density variance, Avg[(XB-XAvg(B))2], where XB is the density of ECs at a given lattice site and XAvg(B) is the average density of ECs taken over the entire simulated space at t = 96 hr (Figure 3C). Here, values close to 0 suggest no distinct separation of the ECs from the SMCs, while positive values indicate varying degrees of separation. As expected, we see a faint positive streak whose position mirrors the zone of co-emergence from the EC fraction plots (Figure 3A). This positive streak approaches 0 with increasing migration values for both ECs and SMCs and aligns with the simulated micropattern mixing observed under high migration values (Figure 3B.iv). Interestingly, micropatterns at lower migration rates, and where cell type mixing is less, were observed to have greater cell-type separation and therefore were more distinct (Figure 3C).
Finally, to quantitatively compare the simulated micropattern predictions with experimentally observed micropatterns, we measured the distributions of EC cluster diameters (see Methods). We found that the majority of the simulated EC cluster diameters, within the physiological range of migration, had effective diameters between 300–800 μm (Figures 3D and S4D). Here, the boundary for the physiological range, JB = 0–0.013 and JC = 0–0.078, outlined by white dashed lines, is set by the maximum reported migration values for ECs [33] and SMCs [42] (Figure 3D). The experimentally observed EC cluster diameters were found to be between 340 ± 110 μm (Figure S1) indicating that our simulations do account for the experimentally observed patterning within the explored physiological range. Last, the simulated EC cluster distribution diameters within the corresponding zone of co-emergence increased as migration increased for both cell types. At very high migration rates, well outside the experimental range, clusters with effective diameters exceeding 800 μm emerged (Figures 3D and S4D), though patterning became less distinct (as seen in Figure 3Biv).
Next, we performed a parameter sweep of proliferation rates in the absence of any sensing mechanisms (Figure 4A). Proliferation rates were varied for ECs (B cells, δB) and SMCs (C cells, δC), while the proliferation rate for VPCs (A cells, δA) and differentiation and migration rates for all cells were fixed at their experimentally relevant values (see Table 1). Reported doubling times include 19 hrs for mouse ECs [43] and 22 hrs for rat vascular SMCs [44]. The parameter sweeps therefore explored physiologically reasonable proliferation rates ranging from no cell division (δθ = 0) to divisions occurring as fast as every 10 hrs (δθ = 0.1). This region in parameter space is demarcated by white dashed lines in Figure 4A. Similar to the phase behavior with varying cell migration, the parameter sweeps here revealed that roughly equal rates of proliferation for ECs, δB, and SMCs, δC, lead to the formation of a zone of co-emergence within which micropatterning develops (Figure 4B i-iii). At the center of the zone, where δB and δC are both 0.05, we found highly distinct micropatterning emerge (Figure 3B ii). As expected, co-culture patterning largely developed for proliferation rates within the confines of the zone of co-emergence while parameter values outside the zone led to the faster growing cell dominating (Figures 4B iv-vi and S4E).
We then looked at the EC density variance plots which shows regions where distinct micropatterns of ECs emerged as a function of proliferation rates (Figure 4C). It was observed that varying the proliferation rates allowed for a larger span of conditions where cell-type separation, defined as having a finite EC cluster density, and thus micropatterning would emerge. Similar to the case with migration, we observe enhanced separation within the zone of co-emergence (Figure 4B i-iii). Interestingly, we observed that some separation persisted beyond the zone of co-emergence, where conditions favored EC proliferation rates over SMC rates (Figure 4C). In these instances, it was observed that the emerging EC clusters were predominately islands of various sizes, and while SMCs were present, they did not in fact surround EC clusters. Additionally, a region of high phase separation was observed at low SMC proliferation rates, δC, and EC proliferation rates of 33 hr, δB = 0.03. This region, which again is outside the zone of co-emergence, consisted largely of only ECs clusters (Figure 4B vi, Figure 4C yellow region).
Last, we measured the distribution of EC cluster diameters and found that they consisted of clusters whose effective diameters were roughly between 500–800 μm, within the physiological range (white dotted lines, Figure 4D), with a broad distribution with standard deviation of 300–800 μm (Figure S4H). Finally, as proliferation rates increased for both ECs and SMCs, a narrow region of parameter space, with EC clusters between 500 and 800 μm, emerged and mirrored the position of the zone of co-emergence (Figure 4D).
Next, we explored the role that differentiation plays in the emergence of micropatterns in the absence of any sensing or signaling. The parameter sweeps were conducted by varying the differentiation rates for both ECs (B cells, αB) and SMCs (C cells, αC) ranging between no differentiation (αθ = 0) and a differentiation time of 10hrs (αθ = 0.1), while holding the migration and proliferation rates constant and at experimentally relevant values (See Table 1, Figure 5A). Narrowing down the physiologically relevant parameter space for differentiation proved to be a difficult task given the variance under different conditions (e.g., chemical, mechanical, and/or contact mediated). Therefore, we define the physiologically relevant differentiation rate range as that between the first/early marker expression to the time point where the given marker expression peaks, indicating a mature phenotype. For mouse ECs it has been reported that VE-cadherin, a known EC marker, expression can be observed as early as 24 hrs post induction [45], corresponding an upper limit of αB = 0.042. Additionally in one of our previous study, we showed that VE-cadherin peaks in differentiating mouse cells around day 14 post induction [25], corresponding a lower limit of αB = 0.003. For mouse smooth muscle cells, early marker expression of α−smooth muscle actin (α-SMA) has been reported at 36 hrs post induction [46], thus establishing an upper limit of αC = 0.028. Furthermore, peak α-SMA expression is reported to arise on day 15 [60], corresponding to a lower limit of αC = 0.0028.
Upon analysis of the EC fraction parameter sweep, we noticed the absence of a sharply defined zone of co-emergence (Figure 5A), defined as EC fractions roughly between 0.3 to 0.7, and instead saw a large fanned out area of parameter space exhibiting micropattern formation (Figure 5B ⅰ-ⅲ). As expected, when both differentiation rates, αB and αC, are zero only VPC (A cells) are present (Figure 5B ⅳ). Additionally, a loss of patterning only occurred when the differentiation rate of one cell type exceeded the other cell type by more than a factor of 2–3 (Figure 5B ⅴ-ⅵ). This is also reflected in the asymmetry plots (Figure S4I) and in the EC density variance plots (Figure 5C), where an increasing gradient of EC density variance values is observed rather than the localized regions of cell-type separation observed in the migration and proliferation EC density plots. This gradient highlights the lack of impact that differentiation rates ultimately have on pattern formation as long as the differentiation rates of EC/SMC cell types are within a factor of 2–3 of each other. Last, the EC cluster diameter distribution was also observed to be fairly constant with cluster diameters ranging from 350 to 800 μm for the total explored parameter space (Figure 5D). A similarly broad range was recorded for the standard deviation, 10 to 500 μm (Figure S4L). To further explore the lower limits of the differentiation rates we repeated our simulations focusing on these longer differentiation periods. Here the range of interest was between 0.001 to 0.01, corresponding to differentiation rates of 1000 hrs [~42 days] to 100 hrs [~4 day]. Unsurprisingly, we did not see any defining differences at these lower bounds (Figure 5E–G).
So far, the results suggest that it is the post-differentiation migration and well-balanced proliferation of the cells that lead to cell patterning. A defining feature of this micropatterning that emerges within the zone of co-emergence is the size of the EC clusters. To understand the effect that migration and proliferation have on a cluster's growth and size, we focused on the growth dynamics of a single EC cluster. We initialized the simulation with 10 cells at the center of the simulation site. As the ECs migrate and proliferate into neighboring lattice sites, the extent of their spreading can be calculated and visualized over time. To do this, we constructed a one-dimensional visualization of the EC spread over time by plotting the EC population along the x-axis as a function of time (Figure 6A). For physiological rates, i.e., when JB = 0.0079 and δB = 0.025, we observed a rapid exodus from the center lattice site into the vacant neighboring lattice sites over time. Additionally, we observed two different regions emerge in the growing EC cluster. First, a leading diffusive front consisting of fractional EC densities, characterized by densities between 0.1 and 1 cells per lattice site. Second, a growing/expanding inner core region in the center, characterized by the presence of 2 or more ECs per lattice site (Figure 6B inset).
To quantify the size of the cluster, we then calculated the radius of gyration of the cluster in the x-direction, defined as σ2 = ∑(xi2*PEC(xi)) where xi is the lattice position and PEC(xi) is the normalized EC population density along the x-axis such that PEC(xi) is equal to XB(xi) / ∑(XB(xi)). The radius of gyration was then plotted over time (Figure 6B–E) followed by a linear and/or quadratic fit, representing diffusive and ballistic regimes respectively, with corresponding R2 values Tables S1–S4). Here, we considered four different cases (ⅰ) physiological migration and proliferation rates (where JB = 0.0079 and dB = 0.025), (ⅱ) high migration (JB = 0.079) with a physiological proliferation rate (dB = 0.025), (ⅲ) physiological migration (JB = 0.0079) with a high proliferation rate (dB = 0.1), and finally (ⅳ) high migration (JB = 0.1) with a high proliferation rate (dB = 0.1).
For case (ⅰ), with both physiological rates, a linear fit was sufficient to describe the time dependence of the radius of gyration (Figure 6B; R2 value of 0.997, Tables S1 and S2). Using the slope of the linear fit, we were able to determine the diffusive rate of the spreading EC cluster to be 281 μm2/hr. Taking the diameter of the cluster to be approximated by 4 σ (which contains 95% of the cells for a Gaussian distribution) allows for direct comparisons with cluster diameters from simulations and experiments that were obtained by binarizing and thresholding images. In this case, we obtain a cluster with diameter 633 μm after 96 hrs. It is to be noted that this diameter is for a single cluster, in the absence of other competing cell types, is much larger than the measured values from experiments and simulations of about 350 μm. This suggests that the leading diffusive front region of the cluster with low cell densities is potentially outcompeted by the surrounding cell type with only the inner core region surviving as a cluster in the competitive environment. Indeed, the inner core region for these physiological parameter values is roughly half the diameter of the full cluster (Figure 6B) yielding a radius of about 320 μm, consistent with measured cluster diameters in the experiment.
For high (10 times greater) migration rates and physiological proliferation rate, we find that a linear fit can explain the radius of gyration time dependence (R2 = 1). As expected, the diffusive constant is roughly 10 times greater than the physiological conditions' at 2931 μm2/hr. The measured cluster diameter at 96 hrs is also much larger at 2108 μm (Figure 6C). Moreover, we observed that the growing cluster does not retain its dense EC inner core, but rather the whole cluster is completely composed of fractional cell densities (lattice sites containing less than one EC). Consequently, while the increase in migration does ultimately produces a larger cluster, the lack of the core as visualized in the cluster by the absence of a bright red area (Figure 6C inset) indicates there may not be any distinct clusters in the competitive environment with multiple cell types. This is consistent with the lack of distinct patterning reported at high migration rates in our simulations (Figure 3B).
Next, we explored physiologically typical migration rates with high (4 times greater) proliferation rates. Here we noticed that the radius of gyration time dependence (Figure 6D) could not be fit by linear function and was better fit by a quadradic (R2 = 0.999) (Tables S3 and S4). From this quadratic fit we could extract the ballistic growth speed of the cluster, which we found to be 15.3 μm/hr. The ballistically growing cluster was observed to display both a smaller leading edge composed of fractional ECs and a larger growing inner core composed of saturated EC lattice sites (bright red area within the cluster consisting of 10 ECs per lattice site Figure 6D inset). This cluster was calculated to have a diameter of 1339 μm (compare to diameter of 633 μm in the physiological case) demonstrating that proliferation alone can drive the growth of the cluster significantly. The presence of a dominant inner core region suggests that high proliferation can also lead to clusters that are distinct in situations where both cell types are present.
Last, for both high migration and high proliferation rates, both linear and quadradic regimes (Figure 6E) were observed, indicating that a transition occurs in the growth dynamics of the cluster, specifically from linear diffusive spread to ballistic growth. We estimated the transition time point as the time at which the R2 value of the linear fit falls below 0.99, which in this case was at 49 hrs (Table S5). A quadratic fit on the remainder of the radius of gyration curve resulted in an R2 value of 0.9996 (Tables S3 and S4). This suggests that high proliferation rates can act as an additional driving force leading to an increasing in cluster size, while migration can synergistically enhance the ballistic growth phase. This is apparent when comparing the cluster size at various time points. At 49 hrs, just after the diffusive growth state, the cluster size is 1647 μm, while at 96 hrs the cluster is at 3862 μm.
In typical cell culture systems, cells sense and often adhere to their neighboring cells via integrin and cadherin binding proteins. These types of cellular interactions have been shown to regulate the cell's migration through contact-mediated inhibition [50]. Therefore, we incorporated these effects into our model by appropriately modifying the migration rates for ECs (B cells) and SMCs (C cells) to allow for specific (Sp) and nonspecific (Nsp) cell adhesions (Eqs (5) and (6)). We examined the three possible distinct combinations of these dependencies: BSp-CSp, BNsp-CSp (same as its inverse), and BNsp-CNsp (Figure 2D). The results show that when both ECs and SMCs display specific adhesions (such as BSp-CSp) the zone of co-emergence, defined as EC fractions between 0.3 and 0.7, is blurred and broadened (Figure 7A) compared to the control case with no contact-inhibition (Figure 3A). At higher migration rates, the contour lines exhibit a slight concave curve due to these cells' migration being significantly slowed by their cell-to-cell adhesions (Figure 7A). Indeed, when one cell type displays specific adhesions and the other displays nonspecific adhesions (BNsp-CSp), the line of co-emergence curves towards the nonspecific adhering cell type (Figure 7B), indicating that the non-specifically adhering cells are more constrained. Last, when both cell types display nonspecific cell adhesions (BNsp-CNsp) they are similarly slowed indicated by the resulting linear contour lines (Figure 7C). These trends are also seen in the corresponding asymmetry plots (Figure S6A–C). The EC density variance (Figure S6G–I) shows fairly distinct micropatterning is present within the zone of co-emergence in all cases except when both cell types have high migration values. The patterning is dramatically more distinct in the zone of co-emergence when one cell type displays specific adhesions and the other displays nonspecific adhesions (Figure S6H). Thus, certain combinations of sensing mechanisms can actually increase the robustness of the patterning.
It is well-known that highly confluent monolayers of many cell types will stop proliferating due to contact inhibition, and a corresponding increase in mechanical constraints [41]. We implement contact-inhibited cell growth in our model by modifying the cell proliferation rates based on the local cell density with either specific or nonspecific sensing, Eqs (7) and (8). The behavior of the EC (B cell) fraction over the explored proliferation rate parameter space was examined for three distinct combinations of specific and nonspecific conditions, BSp-CSp, BNsp-CSp, and BNsp-CNsp (Figure 7D–F), just as explored in cell migration. Adding specific and non-specific sensing to varying proliferation rates was qualitatively similar to those obtained from varying migration, although with stronger affects. When both cells display specific sensing, the zone of co-emergence, defined as the area between 0.3 to 0.7 EC fraction, is again blurred and broadened (Figure 7D) compared to the unconstrained control (Figure 4A). Additionally, when one cell type displays nonspecific sensing and the other displays specific sensing, the line of co-emergence curves towards the nonspecific sensing cell type (Figure 7E). Last, when both cell types display nonspecific sensing and proliferation is similarly slowed down, the contour lines straighten (Figure 7F). These trends can also be seen in the cell asymmetry plots (Figure S6D–F). Like the case with migration, the EC density variance (Figure S6J–L) again shows distinct micropatterning is present within the zone of co-emergence in all cases (except when both cell types have high proliferation values), and particularly so when one cell type displays specific adhesions and the other displays nonspecific adhesions (Figure S6K).
Moreover, we looked at how the patterning depends on cell differentiation rates, α, that are affected by the signaling from neighboring cells, i.e., paracrine signaling. Our model incorporates both same cell-directed differentiation and alternate cell-directed differentiation (see Eq (9) and Figure 2E). This sensing combination and the degree of amplification or suppression is controlled by the magnitude and sign of β (Figure S3). We explored different combinations of paracrine signaling dependence and found that, regardless of the sign and magnitude of β, the EC fraction, over the explored differentiation rate parameter space, was largely unaffected with no significant difference in the formation of micropatterning (Figure 7G–J). This is consistent with the lack of impact that differentiation rates ultimately have on pattern formation as long as the differentiation rates of EC/SMC cell types are similar. For full results see Supplementary Figure S7 and corresponding asymmetry and variance plots (Figure S8)
To compare the relative influence of migration, proliferation, and differentiation on micropattern formation, we explored the EC (B cell) fraction as a function of physiologically relevant ratios of these rates between the different cell types (Figure 7K–M). We first explored the EC fraction as a function of the relative migration rates, JB/JC, and relative differentiation rates, αB/αC (Figure 7K). For physiologically relevant scenarios with ratios between 0.5 to 2, the EC fraction contours run approximately parallel to the αB/αC axis, thus indicating that the emergence of micropatterns is much more sensitive to the ratios of migration rates. Similar results are obtained when examining the proliferation rate ratio, δB/δC, versus the differentiation rate ratio αB/αC (Figure 7L), with the EC fractions dependent on the proliferation more significantly than differentiation. Then, looking at the relative effects of ratios of migration and ratios of proliferation rates for ECs and SMCs (Figure 7M) reveals comparable sensitivity to relative changes in both migration and proliferation rates (Figure S9). The slopes of the EC fraction contour lines are just around ~0.5 indicating a slightly higher sensitivity to the migration rates as compared to the proliferation rates.
Here, we presented an on-lattice stochastic population-based model that qualitatively and quantitatively reproduces the observed 2D micropatterns that emerge during the co-differentiation of ECs and SMCs from VPCs. Our model enables the spatial visualization of this dynamic system over time and examines the effects that various biological processes on cellular micropattern development. Specifically, our computational model explored the roles that cell migration, proliferation, and differentiation, as well as contact inhibition and adjacent cell paracrine signaling have on pattern formation.
Our main finding is that cell migration and proliferation are the key factors driving 2D micropattern formation. While VPC differentiation into ECs and SMCs is required, it is not a predominant driving force in pattern formation. This finding supports an assembly mechanism for the emergence of micropatterns consistent with literature demonstrating that migration [61] and proliferation [38] can similarly drive pattern formation. These two biological processes are also the two main driving forces behind wound healing [35,62], during which, proper organization and patterning of multi-cellular tissues must be executed flawlessly. Furthermore, while some stem cell populations are present within the skin's stem cell niche [63], their role is mainly to supply differentiated cells that will migrate into the wound and proliferate.
The emergence of distinct micropatterns is achieved over a broad area of parameter space, termed the "zone of co-emergence", where EC fractions are roughly between 0.30–0.70. Coupled with our finding that the inclusion of neighbor cell sensing mechanisms impacts the shape of the zone of co-emergence, this proved to be a strong predictor for the development of micropatterns that may be application to other biological systems with implications in directing in-vitro 3D organ morphogenesis.
While cell proliferation and migration rates are the primary driving forces that enable micropatterning, the specific rates and the presence of contact inhibition predict the degree of micropattern development and the length scales of the patterns. Our work on single cluster growth dynamics suggests that the cluster sizes are set by how far the inner core regions, containing high cell densities, can expand and occupy space before being outcompeted by the other cell type. While contact inhibition mitigation of migration and proliferation rates can become important, fewer differences are observed between homotypic and nonspecific sensing mechanisms at low migration and proliferation rates. As rates increase, we see a broad and symmetrical expansion of conditions enabling micropattern development. This is true when ECs and SMCs are both sensing homotypically and nonspecifically. However, when one cell type is sensing homotypically and the other is sensing nonspecifically, there is a reduction in the range of conditions that enable micropattern development accompanied by a significant increase in the distinctness of patterning. This suggests that sensing can help increase the robustness of patterning.
A limitation of our model is the lack of single cell specificity, such as cell shape and polarity, which are important when modeling directional propulsion and cell-to-cell adhesions [64,65]. This level of single cell modeling, at high resolutions, would be excessive for our population-based questions. However, our modeling framework does allow for the modification of our lattice site dimensions, which serves to increase or decrease the overall resolution of our system, within reason, by directly affecting the number of cells that can inhabit a single lattice site. Our modeling framework could also be applicable to study the development of other cells or tissues like the emergence of keratinocytes during differentiation, wound healing, and normal skin repair [66] or the precise staggered patterning of R8 photoreceptor precursors with accessory cells in the neural epithelium of Drosophila eyes [67].
We have explored a range of migration, proliferation, and differentiation values to assess the impact each has on the emergence of multicellular micropatterns within developing vascular tissue and found that patterning is dominated by migration and proliferation under physiologically relevant conditions for VPCs. Our work strongly suggests that, even in the absence of any specific mechanisms that drive segregation, like chemotaxis, 2D micropatterning can emerge as long as cellular fractions are maintained within 0.30–0.70, with cluster sizes being set by the growth rate of single cluster inner cores. Our results suggest that even though micropatterning can occur in the absence of sensing, the presence of such mechanisms greatly increases the robustness of patterning, which could be critical to fidelity in tissue development in the naturally noisy and heterogeneous in vivo setting.
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
We would like to acknowledge useful conversations and guidance from Katie Copenhagen, David Quint, and Lian Wong. This work was funded largely through the NSF-CREST: Center for Cellular and Bio-molecular Machines at UC Merced (NSF-HRD-1547848 and NSF-HRD-2112675), as well as the NSF-Science and Technology Center (STC) for the Emergent Behavior of Integrated Biological Systems (EBICS) Award # 0939511 with partial support from the NSF Center for Engineering Mechanobiology grant CMMI-154857. C = Computational power was provided by the Multi-Environment Computer for Exploration and Discovery (MERCED) cluster at UC Merced (NSF-ACI-1429783). AG also acknowledges support from NSF-DMS-1616926 and the hospitality of the Aspen Center for Physics, which is supported by National Science Foundation Grant PHY-1607611.
The authors declare there is no conflict of interest.
[1] |
H. C. Ott, T. S. Matthiesen, S. K. Goh, L. D. Black, S. M. Kren, T. I. Netoff, et al., Perfusion-decellularized matrix: using nature's platform to engineer a bioartificial heart, Nat. Med., 14 (2008), 213–221. https://doi.org/10.1038/nm1684 doi: 10.1038/nm1684
![]() |
[2] |
J. M. Wainwright, C. A. Czajka, U. B. Patel, D. O. Freytes, K. Tobita, T. W. Gilbert, et al., Preparation of cardiac extracellular matrix from an intact porcine heart, Tissue Eng. Part C Methods, 16 (2010), 525–532. https://doi.org/10.1089/ten.tec.2009.0392 doi: 10.1089/ten.tec.2009.0392
![]() |
[3] |
A. B. Daly, J. M. Wallis, Z. D. Borg, R. W. Bonvillain, B. Deng, B. A. Ballif, et al., Initial binding and recellularization of decellularized mouse lung scaffolds with bone marrow-derived mesenchymal stromal cells, Tissue Eng. Part A, 18 (2012), 1–16. https://doi.org/10.1089/ten.tea.2011.0301 doi: 10.1089/ten.tea.2011.0301
![]() |
[4] |
J. Cortiella, J. Niles, A. Cantu, A. Brettler, A. Pham, G. Vargas, et al., Influence of acellular natural lung matrix on murine embryonic stem cell differentiation and tissue formation, Tissue Eng. Part A, 16 (2010), 2565–2580. https://doi.org/10.1089/ten.tea.2009.0730 doi: 10.1089/ten.tea.2009.0730
![]() |
[5] |
H. C. Ott, B. Clippinger, C. Conrad, C. Schuetz, I. Pomerantseva, L. Ikonomou, et al., Regeneration and orthotopic transplantation of a bioartificial lung, Nat. Med., 16 (2010), 927–933. https://doi.org/10.1038/nm.2193 doi: 10.1038/nm.2193
![]() |
[6] |
K. H. Nakayama, C. A. Batchelder, C. I. Lee, A. F. Tarantal, Decellularized rhesus monkey kidney as a three-dimensional scaffold for renal tissue engineering, Tissue Eng. Part A, 16 (2010), 2207–2216. https://doi.org/10.1089/ten.tea.2009.0602 doi: 10.1089/ten.tea.2009.0602
![]() |
[7] |
G. Orlando, C. Booth, Z. Wang, G. Totonelli, C. L. Ross, E. Moran, et al., Discarded human kidneys as a source of ECM scaffold for kidney regeneration technologies., Biomaterials, 34 (2013), 5915–25. https://doi.org/10.1016/j.biomaterials.2013.04.033 doi: 10.1016/j.biomaterials.2013.04.033
![]() |
[8] |
B. E. Uygun, A. Soto-Gutierrez, H. Yagi, M. L. Izamis, M. A. Guzzardi, C. Shulman, et al., Organ reengineering through development of a transplantable recellularized liver graft using decellularized liver matrix, Nat. Med., 16 (2010), 814–820. https://doi.org/10.1038/nm.2170 doi: 10.1038/nm.2170
![]() |
[9] |
A. Soto-Gutierrez, L. Zhang, C. Medberry, K. Fukumitsu, D. Faulk, H. Jiang, et al., A whole-organ regenerative medicine approach for liver replacement, Tissue Eng. Part C Methods, 17 (2011), 677–686. https://doi.org/10.1089/ten.tec.2010.0698 doi: 10.1089/ten.tec.2010.0698
![]() |
[10] |
Y. H. Hsu, M. L. Moya, C. C. W. Hughes, S. C. George, A. P. Lee, A microfluidic platform for generating large-scale nearly identical human microphysiological vascularized tissue arrays, Lab Chip, 13 (2013), 2990. https://doi.org/10.1039/c3lc50424g doi: 10.1039/c3lc50424g
![]() |
[11] |
K. Sakaguchi, T. Shimizu, S. Horaguchi, H. Sekine, M. Yamato, M. Umezu, et al., In vitro engineering of vascularized tissue surrogates, Sci. Rep., 3 (2013), 1316. https://doi.org/10.1038/srep01316 doi: 10.1038/srep01316
![]() |
[12] |
J. Zhang, L. Chu, Z. Hou, M. P. Schwartz, T. Hacker, V. Vickerman, et al., Functional characterization of human pluripotent stem cell-derived arterial endothelial cells, Proc. Natl. Acad. Sci., (2017), 201702295. https://doi.org/10.1073/pnas.1702295114 doi: 10.1073/pnas.1702295114
![]() |
[13] |
J. A. Whisler, M. B. Chen, R. D. Kamm, Control of perfusable microvascular network morphology using a multiculture microfluidic system, Tissue Eng. Part C Methods, 20 (2014), 543–552. https://doi.org/10.1089/ten.tec.2013.0370 doi: 10.1089/ten.tec.2013.0370
![]() |
[14] |
Z. Wan, S. Zhang, A. X. Zhong, S. E. Shelton, M. Campisi, S. K. Sundararaman, et al., A robust vasculogenic microfluidic model using human immortalized endothelial cells and Thy1 positive fibroblasts, Biomaterials, 276 (2021), 121032. https://doi.org/10.1016/j.biomaterials.2021.121032 doi: 10.1016/j.biomaterials.2021.121032
![]() |
[15] |
K. Yamamoto, K. Tanimura, M. Watanabe, H. Sano, H. Uwamori, Y. Mabuchi, et al., Construction of continuous capillary networks stabilized by pericyte-like perivascular cells, Tissue Eng. Part A, 25 (2019), 499–510. https://doi.org/10.1089/ten.tea.2018.0186 doi: 10.1089/ten.tea.2018.0186
![]() |
[16] |
S. Levenberg, J. S. Golub, M. Amit, J. Itskovitz-Eldor, R. Langer, Endothelial cells derived from human embryonic stem cells, Proc. Natl. Acad. Sci., 99 (2002), 4391–4396. https://doi.org/10.1073/pnas.032074999 doi: 10.1073/pnas.032074999
![]() |
[17] |
K. E. McCloskey, D. A. Smith, H. Jo, R. M. Nerem, Embryonic stem cell-derived endothelial cells may lack complete functional maturation in vitro, J. Vasc. Res., 43 (2006), 411–421. https://doi.org/10.1159/000094791 doi: 10.1159/000094791
![]() |
[18] |
Z. Gong, L. E. Niklason, Use of human mesenchymal stem cells as alternative source of smooth muscle cells in vessel engineering, Methods Mol. Biol., 698 (2011), 279–294. https://doi.org/10.1007/978-1-60761-999-4_21 doi: 10.1007/978-1-60761-999-4_21
![]() |
[19] |
G. K. Owens, Regulation of differentiation of vascular smooth muscle cells, Physiol. Rev., 75 (1995), 487–517. https://doi.org/10.1152/physrev.1995.75.3.487 doi: 10.1152/physrev.1995.75.3.487
![]() |
[20] |
E. M. Shen, K. E. McCloskey, Development of mural cells: From in vivo understanding to in vitro recapitulation, Stem Cells Dev., 26 (2017), 1020–1041. https://doi.org/10.1089/scd.2017.0020 doi: 10.1089/scd.2017.0020
![]() |
[21] |
B. Descamps, C. Emanueli, Vascular differentiation from embryonic stem cells: Novel technologies and therapeutic promises, Vascul. Pharmacol., 56 (2012), 267–279. https://doi.org/10.1016/j.vph.2012.03.007 doi: 10.1016/j.vph.2012.03.007
![]() |
[22] |
R. A. Wimmer, A. Leopoldi, M. Aichinger, N. Wick, B. Hantusch, M. Novatchkova, et al., Human blood vessel organoids as a model of diabetic vasculopathy, Nature, 565 (2019), 505–510. https://doi.org/10.1038/s41586-018-0858-8 doi: 10.1038/s41586-018-0858-8
![]() |
[23] |
A. A. Blancas, A. J. Shih, N. E. Lauer, K. E. McCloskey, Endothelial cells from embryonic stem cells in a chemically defined medium, Stem Cells Dev., 20 (2011), 2153–2161. https://doi.org/10.1089/scd.2010.0432 doi: 10.1089/scd.2010.0432
![]() |
[24] |
K. McCloskey, D. Glaser, A. Burns, R. Hatano, Y. Fan, M. Medrzycki, Specialized mouse embryonic stem cells for studying vascular development, Stem Cells Cloning Adv. Appl., 7 (2014), 79. https://doi.org/10.2147/SCCAA.S69554 doi: 10.2147/SCCAA.S69554
![]() |
[25] |
D. E. Glaser, W. S. Turner, N. Madfis, L. Wong, J. Zamora, N. White, et al., Multifactorial optimizations for directing endothelial fate from stem cells, PLoS One, 11 (2016), e0166663. https://doi.org/10.1371/journal.pone.0166663 doi: 10.1371/journal.pone.0166663
![]() |
[26] |
B. Jahan, K. E. McCloskey, Differentiation and expansion of endothelial cells requires pre-optimization of KDR+ expression kinetics, Stem Cell Res., 42 (2020), 101685. https://doi.org/10.1016/j.scr.2019.101685 doi: 10.1016/j.scr.2019.101685
![]() |
[27] |
L. Wong, A. Kumar, B. Gabela-Zuniga, J. Chua, G. Singh, C. L. Happe, et al., Substrate stiffness directs diverging vascular fates, Acta Biomater., 96 (2019), 321–329. https://doi.org/10.1016/j.actbio.2019.07.030 doi: 10.1016/j.actbio.2019.07.030
![]() |
[28] |
J. M. Osborne, A. G. Fletcher, J. M. Pitt-Francis, P. K. Maini, D. J. Gavaghan, Comparing individual-based approaches to modelling the self-organization of multicellular tissues, PLOS Comput. Biol., 13 (2017), e1005387. https://doi.org/10.1371/journal.pcbi.1005387 doi: 10.1371/journal.pcbi.1005387
![]() |
[29] |
D. Viens, A three-dimensional finite element model for the mechanics of cell-cell interactions, J. Biomech. Eng., 129 (2007), 651. https://doi.org/10.1115/1.2768375 doi: 10.1115/1.2768375
![]() |
[30] |
G. W. Brodland, D. A. Clausi, Embryonic tissue morphogenesis modeled by FEM, J. Biomech. Eng., 116 (1994), 146–155. https://doi.org/10.1115/1.2895713 doi: 10.1115/1.2895713
![]() |
[31] |
Q. Smith, E. Stukalin, S. Kusuma, S. Gerecht, S. X. Sun, Stochasticity and spatial interaction govern stem cell differentiation dynamics, Sci. Rep., 5 (2015), 1–10. https://doi.org/10.1038/srep12617 doi: 10.1038/srep12617
![]() |
[32] |
W. C. Lo, C. S. Chou, K. Gokoffski, F. Wan, A. Lander, A. Calof, et al., Feedback regulation in multistage cell lineages, Math. Biosci. Eng., 6 (2009), 59–82. https://doi.org/10.3934/mbe.2009.6.59 doi: 10.3934/mbe.2009.6.59
![]() |
[33] |
A. Szabó, R. Ünnep, E. Méhes, W. O. Twal, W. S. Argraves, Y. Cao, et al., Collective cell motion in endothelial monolayers, Phys. Biol., 7 (2010), 046007. https://doi.org/10.1088/1478-3975/7/4/046007 doi: 10.1088/1478-3975/7/4/046007
![]() |
[34] |
Y. Mao, A. L. Tournier, P. A. Bates, J. E. Gale, N. Tapon, B. J. Thompson, Planar polarization of the atypical myosin Dachs orients cell divisions in Drosophila, Genes Dev., 25 (2011), 131–136. https://doi.org/10.1101/gad.610511 doi: 10.1101/gad.610511
![]() |
[35] |
D. C. Walker, G. Hill, S. M. Wood, R. H. Smallwood, J. Southgate, Agent-based computational modeling of wounded epithelial cell monolayers, IEEE Trans. Nanobiosci., 3 (2004), 153–163. https://doi.org/10.1109/TNB.2004.833680 doi: 10.1109/TNB.2004.833680
![]() |
[36] |
E. L. Bearer, J. S. Lowengrub, H. B. Frieboes, Y. L. Chuang, F. Jin, S. M. Wise, et al., Multiparameter computational modeling of ttumor invasion, Cancer Res., 69 (2009), 4493–4501. https://doi.org/10.1158/0008-5472.CAN-08-3834 doi: 10.1158/0008-5472.CAN-08-3834
![]() |
[37] |
D. Drasdo, S. Höhme, A single-cell-based model of tumor growth in vitro : monolayers and spheroids, Phys. Biol., 2 (2005), 133–147. https://doi.org/10.1088/1478-3975/2/3/001 doi: 10.1088/1478-3975/2/3/001
![]() |
[38] |
C. S. Chou, W. C. Lo, K. K. Gokoffski, Y. T. Zhang, F. Y. M. Wan, A. D. Lander, et al., Spatial dynamics of multistage cell lineages in tissue stratification, Biophys. J., 99 (2010), 3145–3154. https://doi.org/10.1016/j.bpj.2010.09.034 doi: 10.1016/j.bpj.2010.09.034
![]() |
[39] |
D. J. Kelly, P. J. Prendergast, Mechano-regulation of stem cell differentiation and tissue regeneration in osteochondral defects, J. Biomech., 38 (2005), 1413–1422. https://doi.org/10.1016/j.jbiomech.2004.06.026 doi: 10.1016/j.jbiomech.2004.06.026
![]() |
[40] |
H. Du, Y. Wang, D. Haensel, B. Lee, X. Dai, Q. Nie, Multiscale modeling of layer formation in epidermis, PLOS Comput. Biol., 14 (2018), e1006006. https://doi.org/10.1371/journal.pcbi.1006006 doi: 10.1371/journal.pcbi.1006006
![]() |
[41] |
A. Atala, Re: Collective and single cell behavior in epithelial contact inhibition, J. Urol., 188 (2012), 1396–1397. https://doi.org/10.1016/j.juro.2012.06.073 doi: 10.1016/j.juro.2012.06.073
![]() |
[42] |
B. Schreier, G. Schwerdt, C. Heise, D. Bethmann, S. Rabe, S. Mildenberger, et al., Substance-specific importance of EGFR for vascular smooth muscle cells motility in primary culture, Biochim. Biophys. Acta Mol. Cell Res., 1863 (2016), 1519–1533. https://doi.org/10.1016/j.bbamcr.2016.03.017 doi: 10.1016/j.bbamcr.2016.03.017
![]() |
[43] |
J. Walter-Yohrling, S. Morgenbesser, C. Rouleau, R. Bagley, M. Callahan, W. Weber, et al., Murine endothelial cell lines as models of tumor endothelial cells, Clin. Cancer Res., 10 (2004), 2179–2189. https://doi.org/10.1158/1078-0432.CCR-03-1013 doi: 10.1158/1078-0432.CCR-03-1013
![]() |
[44] |
N. Endlich, K. Endlich, N. Taesch, J. J. Helwig, Culture of vascular smooth muscle cells from small arteries of the rat kidney, Kidney Int., 57 (2000), 2468–2475. https://doi.org/10.1046/j.1523-1755.2000.00105.x doi: 10.1046/j.1523-1755.2000.00105.x
![]() |
[45] |
S. I. Nishikawa, S. Nishikawa, M. Hirashima, N. Matsuyoshi, H. Kodama, Progressive lineage analysis by cell sorting and culture identifies FLK1+VE-cadherin+ cells at a diverging point of endothelial and hemopoietic lineages, Development, 125 (1998), 1747–1757. https://doi.org/10.1111/ijpo.259 doi: 10.1111/ijpo.259
![]() |
[46] |
J. Yamashita, H. Itoh, M. Hirashima, M. Ogawa, S. Nishikawa, T. Yurugi, et al., Flk1-positive cells derived from embryonic stem cells serve as vascular progenitors, Nature, 408 (2000), 92–96. https://doi.org/10.1038/35040568 doi: 10.1038/35040568
![]() |
[47] |
K. L. Hill, P. Obrtlikova, D. F. Alvarez, J. A. King, S. A. Keirstead, J. R. Allred, et al., Human embryonic stem cell−derived vascular progenitor cells capable of endothelial and smooth muscle cell function, Exp. Hematol., 38 (2010), 246–257. https://doi.org/10.1016/j.exphem.2010.01.001 doi: 10.1016/j.exphem.2010.01.001
![]() |
[48] | G. Van Rossum, F. L. Drake Jr, Python reference manual, Department of Computer Science, CWI, 1995. |
[49] | MATLAB, version 9.1.0.441655 (R2016b), Natick, MA: The MathWorks, Inc. Available from: https://www.mathworks.com |
[50] |
A. Huttenlocher, M. Lakonishok, M. Kinder, S. Wu, T. Truong, K. A. Knudsen, et al., Integrin and cadherin synergy regulates contact inhibition of migration and motile activity, J. Cell Biol., 141 (1998), 515–526. https://doi.org/10.1083/jcb.141.2.515 doi: 10.1083/jcb.141.2.515
![]() |
[51] |
M. Abercrombie, J. E. M. Heaysman, Observations on the social behaviour of cells in tissue culture, Exp. Cell Res., 6 (1954), 293–306. https://doi.org/10.1016/0014-4827(54)90176-7 doi: 10.1016/0014-4827(54)90176-7
![]() |
[52] |
M. S. Steinberg, Does differential adhesion govern self‐assembly processes in histogenesis? Equilibrium configurations and the emergence of a hierarchy among populations of embryonic cells, J. Exp. Zool., 173 (1970), 395–433. https://doi.org/10.1002/jez.1401730406 doi: 10.1002/jez.1401730406
![]() |
[53] |
B. Lilly, We have contact: endothelial cell-smooth muscle cell interactions, Physiology, 29 (2014), 234–241. https://doi.org/10.1152/physiol.00047.2013 doi: 10.1152/physiol.00047.2013
![]() |
[54] |
W. Baumgartner, P. Hinterdorfer, W. Ness, A. Raab, D. Vestweber, H. Schindler, et al., Cadherin interaction probed by atomic force microscopy, Proc. Natl. Acad. Sci., 97 (2000), 4005–4010. https://doi.org/10.1073/pnas.070052697 doi: 10.1073/pnas.070052697
![]() |
[55] |
E. Moiseeva, Adhesion receptors of vascular smooth muscle cells and their functions, Cardiovasc. Res., 52 (2001), 372–386. https://doi.org/10.1016/S0008-6363(01)00399-6 doi: 10.1016/S0008-6363(01)00399-6
![]() |
[56] |
E. Perret, A. Leung, H. Feracci, E. Evans, Trans-bonded pairs of E-cadherin exhibit a remarkable hierarchy of mechanical strengths, Proc. Natl. Acad. Sci., 101 (2004), 16472–16477. https://doi.org/10.1073/pnas.0402085101 doi: 10.1073/pnas.0402085101
![]() |
[57] |
M. Noseda, L. Chang, G. McLean, J. E. Grim, B. E. Clurman, L. L. Smith, et al., Notch activation induces endothelial cell cycle arrest and participates in contact inhibition: role of p21Cip1 repression, Mol. Cell. Biol., 24 (2004), 8813–8822. https://doi.org/10.1128/MCB.24.20.8813-8822.2004 doi: 10.1128/MCB.24.20.8813-8822.2004
![]() |
[58] |
J. P. Sasine, K. T. Yeo, J. P. Chute, Concise review: Paracrine functions of vascular niche cells in regulating hematopoietic stem cell fate, Stem Cells Transl. Med., 6 (2017), 482–489. https://doi.org/https://doi.org/10.5966/sctm.2016-0254 doi: 10.5966/sctm.2016-0254
![]() |
[59] |
H. L. Kirschenlohr, J. C. Metcalfe, P. L. Weissberg, D. J. Grainger, Adult human aortic smooth muscle cells in culture produce active TGF-beta, Am. J. Physiol., 265 (1993), C571–C576. https://doi.org/10.1152/ajpcell.1993.265.2.C571 doi: 10.1152/ajpcell.1993.265.2.C571
![]() |
[60] |
H. Huang, X. Zhao, L. Chen, C. Xu, X. Yao, Y. Lu, et al., Differentiation of human embryonic stem cells into smooth muscle cells in adherent monolayer culture, Biochem. Biophys. Res. Commun., 351 (2006), 321–327. https://doi.org/10.1016/j.bbrc.2006.09.171 doi: 10.1016/j.bbrc.2006.09.171
![]() |
[61] |
E. G. Rens, M. T. Zeegers, I. Rabbers, A. Szabó, R. M. H. Merks, Autocrine inhibition of cell motility can drive epithelial branching morphogenesis in the absence of growth, Philos. Trans. R. Soc. B Biol. Sci., 375 (2020), 20190386. https://doi.org/10.1098/rstb.2019.0386 doi: 10.1098/rstb.2019.0386
![]() |
[62] |
S. H. Yoon, Y. K. Kim, E. D. Han, Y. H. Seo, B. H. Kim, M. R. K. Mofrad, Passive control of cell locomotion using micropatterns: the effect of micropattern geometry on the migratory behavior of adherent cells, Lab Chip, 12 (2012), 2391. https://doi.org/10.1039/c2lc40084g doi: 10.1039/c2lc40084g
![]() |
[63] |
A. Spradling, D. Drummond-Barbosa, T. Kai, Stem cells find their niche, Nature, 414 (2001), 98–104. https://doi.org/10.1038/35102160 doi: 10.1038/35102160
![]() |
[64] |
P. J. Albert, U. S. Schwarz, Modeling cell shape and dynamics on micropatterns, Cell Adh. Migr., 10 (2016), 516–528. https://doi.org/10.1080/19336918.2016.1148864 doi: 10.1080/19336918.2016.1148864
![]() |
[65] |
X. Jiang, D. A. Bruzewicz, A. P. Wong, M. Piel, G. M. Whitesides, Directing cell migration with asymmetric micropatterns, Proc. Natl. Acad. Sci., 102 (2005), 975–978. https://doi.org/10.1073/pnas.0408954102 doi: 10.1073/pnas.0408954102
![]() |
[66] |
N. Ojeh, I. Pastar, M. Tomic-Canic, O. Stojadinovic, Stem cells in skin regeneration, wound healing, and their clinical applications, Int. J. Mol. Sci., 16 (2015), 25476–25501. https://doi.org/10.3390/ijms161025476 doi: 10.3390/ijms161025476
![]() |
[67] |
R. I. Johnson, Hexagonal patterning of the Drosophila eye, Dev. Biol., 478 (2021), 173–182. https://doi.org/10.1016/j.ydbio.2021.07.004 doi: 10.1016/j.ydbio.2021.07.004
![]() |
![]() |
![]() |
Variable | Denotation | Experimentally obtained values | Corresponding simulation value(s) | Parameter sweep | References |
Diffusion constant | Jθ | 14 μm2/hr | 0.0079 | 0–0.1; 0.001 step size | Supplemental video V1 (37–38) |
Proliferation rate | δθ | ~40 hr | 0.025 | 0–0.1; 0.001 step size | (39–41) |
Differentiation rate | αθ | ~62.5 hr | 0.016 | 0–0.1; 0.001 step size | (19, 20) |
Stochastic noise | Sθ | A loss or gain of up to 1 cell per lattice site is incorporated every 10 hrs | ± 0.01 | N/A | NA |
Paracrine signal strength | β | amplification of differentiation rates by 0.5–1.5 | ± 0.5, ± 0.016, 0 | N/A | NA |
Variable | Denotation | Experimentally obtained values | Corresponding simulation value(s) | Parameter sweep | References |
Diffusion constant | Jθ | 14 μm2/hr | 0.0079 | 0–0.1; 0.001 step size | Supplemental video V1 (37–38) |
Proliferation rate | δθ | ~40 hr | 0.025 | 0–0.1; 0.001 step size | (39–41) |
Differentiation rate | αθ | ~62.5 hr | 0.016 | 0–0.1; 0.001 step size | (19, 20) |
Stochastic noise | Sθ | A loss or gain of up to 1 cell per lattice site is incorporated every 10 hrs | ± 0.01 | N/A | NA |
Paracrine signal strength | β | amplification of differentiation rates by 0.5–1.5 | ± 0.5, ± 0.016, 0 | N/A | NA |