Dynamic intercellular transport modulates the spatial patterning of differentiation during early neural commitment

We developed a multicellular model that describes cell fate transitions coupled with intercellular communication that is regulated by the dynamics of connexins as a function of the cell cycle. Although previous computational models of intercellular communication have included the biophysics of gap junction transport, they have primarily focused on small molecule diffusion through individual channels48,49. One notable multicellular model that focused on ion transfer relied upon membrane potential to regulate intercellular connectivity50. While capable of examining the development of intracellular gradients at a multicellular level, the bioelectric model is constrained to a static system. A primary advantage of our selected framework is its ability to investigate the impact of dynamic intercellular topologies on a multicellular system that includes growth, division, and differentiation. The dynamic network facilitates the study of evolving spatial organization and the resultant effect on the formation of intracellular gradients.

Interestingly, a recent paper proposed an alternative, hidden Markov model framework that explored the role of stochasticity during early neural differentiation, but without consideration of spatial organization and cell-cell communication51. While the mechanistic framework is different from the work presented here, they share several similarities. In particular, both models consider differentiation to be an autonomous process with discrete state changes that are probabilistic. The probability of a state transition in our model, however, is dependent on the accumulation of an intracellular metabolite within a cell rather than the purely stochastic progression through microstates. The description of microstates typically considers gene and protein expression patterns, but in the present work the metabolic profile of a cell also impacts susceptibility to differentiation signals52,53. The connection between microstates and the intracellular accumulation of metabolites is further highlighted by the temporal delay effect of both AC- and GJ-inhibition on the trajectory of differentiation (i.e. slower microstate progression). The cellular communication integrated within our model via gap junction-mediated transport leverages both stochastic state transitions (consistent with Stumpf et al.) and the higher-order collective behavior of the multicellular system, suggested by our observations of pattern propagation within pluripotent colonies. Therefore, our model both complements this prior modeling study while also advancing the description of spatial organization effects on the emergence and propagation of micro- and macrostates during early neural differentiation.

The superposition in latent variable space between experimentally derived results and agent-based simulations is a powerful feature of the analysis pipeline presented here. The extraction of network metrics is particularly necessary for the accurate quantification of systems with variable and asymmetrical structures, such as during unconstrained monolayer growth (Fig. 4c). Since each colony structure can have a unique configuration, the spatial characteristics used to distinguish pattern classes need to be independent from the underlying morphology. Network metrics accomplish this by evaluating relationships between neighboring nodes/cells and colony-wide. Furthermore, the transformation of each multivariate set of metrics to latent space enables visual and numerical comparison between individual colonies. The training set of artificially generated colony patterns provides the ability to directly map regions of latent space to specific spatial organizations (e.g. the “random” cloud versus a highly organized “outside-in” or “inside-out” cloud). Accordingly, transformed sets of colony metrics from both experimental and computational sources can be delineated according to their proximity to each mapped pattern class in latent space. Computational models therefore can be designed to replicate morphogenic trajectories in a quantifiable manner and elucidate the factors that influence specific pattern formations.

In summary, this work demonstrates that intercellular heterogeneity can be generated by dynamic fluctuations in GJ permeability via the cell cycle. We found that the spatial trajectory of early neural differentiation is remarkably robust to perturbation through small molecule inhibition, but that the temporal progression of differentiation can be modulated. Ultimately, identifying the effect of molecular regulators of differentiation at a multicellular scale is an initial step towards replicating morphogenic events in vitro and successfully directing the emergence of desired engineered multicellular structures.

The murine embryonic stem cell line (D3) was cultured at 37 °C in monolayer on 100-mm tissue culture plates coated with 0.1% gelatin (Millipore EmbryoMax) in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 15% fetal bovine serum (FBS) (Atlanta Biologicals, Atlanta, GA), 2 mM L-glutamine (Lonza), 100 U mL−1 penicillin, 100 μg mL−1 streptomycin, and 0.25 μg mL−1 amphotericin (MP Biomedicals), 1x MEM nonessential amino acid solution (Corning), 0.1 mM 2-mercaptoethanol (Sigma-Aldrich), and 103 mL−1 leukemia inhibitory factor (LIF) (EMD Millipore). Cells were passaged every 2–3 days, using 0.05% trypsin (Corning) to dissociate cells, centrifuged at 200 rcf for 5 min, and plated at a density of 20000 cells cm−2.

A MycoAlert™ PLUS Mycoplasma Detection Kit (Lonza LT07-703) was used to verify that there was no mycoplasma contamination in our cultures.

Cells were plated at a density of 10000 cells cm−2on Ibidi μ-slides coated with a 0.1% gelatin, 0.05% fibronectin solution. Cells were allowed to grow under regular cell culture conditions described earlier for 24 h before being transitioned to N2B27 medium containing 1 μM retinoic acid (Sigma-Aldrich). The N2B27 medium was composed of 50% DMEM/F12 (Thermo Fisher Scientific) and 50% Neurobasal Media (Thermo Fisher Scientific), supplemented with 0.5% N-2(100×) (Thermo Fisher Scientific) and 1% B-27(50×) (Thermo Fisher Scientific). Every 24 h, the medium was replaced with fresh N2B27 medium containing 1 μM retinoic acid.

For staining of Oct4 and Sox2, cells were fixed in 4% paraformaldehyde (Thermo Fisher Scientific) for 5 min and then permeabilized with 0.1% Triton X100 for 15 min. Samples were blocked with 2% normal donkey serum (NDS) (Sigma-Aldrich) for 1 hour and then incubated overnight at 4 °C with the following primary antibodies in 2% NDS: goat polyclonal Oct-3/4 (Santa Cruz, 1:200) and rabbit polyclonal Sox-2 (Thermo Fisher Scientific, 1:200). After washing, samples were incubated in a 2% NDS secondary antibody solution of donkey anti-goat Alexa Fluor 568 (Thermo Fisher Scientific, 1:200) and donkey anti-rabbit Alexa Fluor 647 (Thermo Fisher Scientific, 1:200) for 45 min before counterstaining with Hoechst 33342 (1:1000) in DI water. All images were collected on a PerkinElmer UltraVIEW VoX spinning-disk confocal microscope with a sCMOS camera at 20X magnification.

For Cx43 and Oct4 staining, the same procedure was utilized but with the following changes: The duration of fixation was increased to 15 min, the primary antibodies were goat polyclonal Oct-3/4 (Santa Cruz, 1:200) and rabbit polyclonal Cx43 (Sigma, 1:400), and the images were collected at 60x magnification. Catalogue numbers for antibodies are in Table 1.Table 1List of antibodies and their respective catalog numbersTargetCompanyCatalog #HostOct-3/4 (N-19)Santa Cruzsc-8628GoatSOX2InvitrogenPA1-16968RabbitCx43Sigma-AldrichC6219Rabbit

To form colonies, mouse (m)ESCs were plated at a density of 10,000 cells cm−2 on Ibidi μ-slides, coated with a 0.1% gelatin, 0.05% fibronectin solution, for 24 h. A 1 mM calcein-AM (Thermo Fisher Scientific) stock solution in DMSO was diluted into 100 μL DMEM media and vortexed. The diluted calcein-AM solution was added to each well of the Ibidi μ-slide to a final concentration of 1 μM for 40 min at 37 °C. After 40 min, Alexa Fluor 647–conjugated wheat germ agglutinin (WGA) (Thermo Fisher Scientific) was added to each well to a final concentration of 10 μg mL−1, for 5 min at 37 °C. Each well was rinsed once before phenol red–free DMEM (Corning) medium was added for imaging. Images were acquired with a Zeiss NLO 710 confocal microscope, using built in software for photobleaching individual cells. Bleaching was initiated after acquiring three images, then maximum laser intensity was pulsed for 30 iterations before measuring fluorescence intensity in the photobleached cell every 3.8 sec. For identifying mitotic cells, Hoechst 33342 was added with the WGA for 5 min at a 1:1000 dilution. When quantifying nocodazole-treated cells, each well was treated with 100 ng mL−1 of nocodazole (Sigma-Aldrich) for 4 h before replacing with new DMEM medium containing 1 μM calcein-AM and incubating for 40 min at 37 °C, followed by 5 min with WGA, then changing to phenol red–free DMEM for imaging.

From each GAP-FRAP experiment, the measured fluorescence intensity value in the photobleached area was compiled as a function of time, setting t = 0 as the time point with the minimum intensity value after bleaching. The intensity values were normalized to the maximum intensity after bleaching and fitted to an exponential function of the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I = 1 - {\mathrm{exp}}\left( { - \frac{t}{{R_{\mathrm{c}}}}} \right)$$\end{document}I=1-exp-tRc where Rc is the recovery constant. The recovery constant was then normalized to the maximum fluorescence intensity before bleaching. In a few cases, it was also possible to quantify the recovery lost from neighboring donor cells (Supplementary Figures 1, 4). However, the context of the intercellular transport is different for the neighboring cells compared to the photobleached cell. Specifically, the photobleached cell only experiences a net influx of unbleached calcein, whereas the neighboring cells have both a loss of calcein to the photobleached cell and an influx from secondary neighbors. As such, measurements from secondary cells are not included in any of the distributions for characterizing transport rates in Fig. 3.

For quantifying differentiation progression, cells were grown for 24 h on 100 mm tissue culture plates coated with 0.1% gelatin. Each plate was rinsed with PBS then differentiated using our previously stated differentiation protocol, with samples taken at 0, 24, 48, and 72 h from the start of the protocol. Each sample was dissociated using 0.05% trypsin, centrifuged for 5 min at 200 rcf, and fixed using 4% paraformaldehyde for 10 min. Cells were permeabilized by centrifuging at 200 rcf with 5% Triton X100 for 5 min. To prevent non-specific binding, cells were incubated for an hour in 10% NDS before being treated with primary antibodies Oct-3/4 (Santa Cruz, 1:200) and Sox-2 (Thermo Fisher Scientific, 1:200) in 10% NDS for 1 h. Cells were washed 3 times then incubated for 30 min in 10% NDS containing the following secondary antibodies: Donkey anti-goat Alexa Fluor 488 (Thermo Fisher Scientific, 1:200) and Donkey anti-rabbit Alexa Fluor 647 (Thermo Fisher Scientific, 1:200).

For cell cycle analysis, the same protocol was used but with incubation in Hoechst 33342 (1:1000) for 15 min instead of primary and secondary antibody staining. All flow cytometry measurements were collected on a BD LSR Fortessa flow cytometer.

An agent-based model was developed for simulation of pluripotent colonies by refining a python-based custom platform first described in White et al33. Here, each cellular agent with a set of properties of cell radius, spatial coordinates, division time, cell type, and intracellular concentration. Cell radius was set at 6.5 µm based on previous data collected in our lab for D3 cells, with cells modeled as rigid bodies. Spatial coordinates were initially set using a digitized set of coordinates from an experimental colony, assuming 100% pluripotent cells. Division time is randomly distributed between 0 and 18 h for each cell and is defined as the time since the cell last divided, with pluripotent cells dividing every 18 h and differentiated cells dividing every 51 h. As the colony grows it is restricted to a 2D-plane, with a collision detection algorithm optimizing cellular movement to prevent overlap of cells in the model. The asymmetrical growth of the colony is accomplished by preferential movement of cells away from the central mass of the colony when resolving overlap during collision detection. The initial intracellular concentration of each cell is randomly selected between 0.6 and 0.65, with a nascent permeability, production rate, and degradation rate specific to each cell type (Pluripotent, Differentiated). Table 2 presents parameter values for each cell type.Table 2Cell type specific parameters used in ABM of intercellular communicationCell typePluripotentDifferentiatedCell radius (r)6.56.5Cell cycle length (CCL)18 h51 hNascent permeability (PMn)0.450.85Production constant (vp)4.89E-66.0E-6 K 50 0.030.04Degradation constant (Kd)1.3E-71.3E-7

The intercellular network is determined at each time step by calculating the intercellular distance between every cell, and if the cell membranes are within 2 microns (2xradius + 2 um) then they are “connected”. Since the resistance to diffusion is significantly higher at the gap junction interface, we consider the intracellular compartment to have instantaneous mixing after every diffusion time step. The flux between each cell follows Fick’s first law, accounting for the individual permeability of the connected cells (function of cell cycle and cell type, Eq. 1–6.

PMmax: Maximum permeability, the maximum flux between any two cells if all gap junctions were open1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{DT}}_{{\mathrm{norm}}} = {\mathrm{Abs}}\left( {{{t}}_{{\mathrm{Div}}} - \left( {\frac{{{\mathrm{CCL}}}}{4}} \right)} \right)$$\end{document}DTnorm=AbstDiv-CCL4

DTnorm is a transform function, shifting the minimum value forward by ¼ of the cell-cycle length (Supplementary Figure 1,5). This shift decreases the steep transition in cell-cycle time that is seen in tdiv when a cell starts a new division cycle (@18 h simulation time), while maintaining individual regions of minimum and maximum values per period.2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{Cell}}{\hskip2pt}{\mathrm{Cycle}}{\hskip2pt}{\mathrm{Effect}}\left( {{\mathrm{CCE}}} \right) = \frac{1}{{1 + \left( {\frac{{{\mathrm{DT}}_{{\mathrm{norm}}}}}{{{\mathrm{CCL}}}}} \right)^2}} + \frac{{0.69t_{{\mathrm{Div}}}^6}}{{{{t}}_{{\mathrm{Div}}}^6 + {\mathrm{CCL}}^6}}$$\end{document}CellCycleEffectCCE=11+DTnormCCL2+0.69tDiv6tDiv6+CCL6

The cell-cycle effect function is defined such that maximum transport would occur during the first few h of a new division cycle when a cell would be in G1-phase, followed by a decline throughout S-phase to reach a minimum transport value near the end of a division cycle when a cell would be undergoing mitosis.3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{PM}}_{{\mathrm{n},{\mathrm {eff}}}} = {\mathrm{PM}}_{{{\max}}} \left( {\mathrm{CCE}}_{{\mathrm{cell 1}}} \ast {\mathrm{PM}}_{{\mathrm{n}},{\mathrm{cell}}1} \right)\left( {\mathrm{CCE}}_{\mathrm{cell 2}} \ast {\mathrm{PM}}_{{\mathrm{n}},{\mathrm{cell 2}}} \right)$$\end{document}PMn,eff=PMmaxCCEcell1*PMn,cell1CCEcell2*PMn,cell2

The effective permeability between any two cells is equal to the maximum permeability adjusted by the product of the nascent permeability of the adjacent cells. The nascent permeability is defined by cell type (differentiated or undifferentiated) and regulated by the Cell cycle effect function.4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{Prod}} = \frac{{{\mathrm{vp}}}}{{1 + {{C}}_{\mathrm{t}}/K_{50}}}$$\end{document}Prod=vp1+Ct∕K505\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{Deg}} = - k_{\mathrm{d}}C_{\mathrm{t}}$$\end{document}Deg=-kdCt6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{Flux}} = {\mathrm{PM}}_{{\mathrm{n}},{\mathrm{eff}}}\Delta C + {\mathrm{Prod}} + {\mathrm{Deg}}$$\end{document}Flux=PMn,effΔC+Prod+Deg

Production, degradation, and flux in and out of a cell are calculated every 3 s, or 1200 times per differentiation time step, with a continuously updating intracellular concentration. For the production function, we use an inhibitory Hill function with an n-value of 1 to represent decreased production when intracellular concentrations are high. A simple degradation function is used with the degradation constant the same for both cell types. Diffusive flux between adjacent cells is calculated as the product of the effective permeability and the concentration gradient between each neighboring cell. The final concentration at the end of every hour is used to determine the probability of differentiation, as per Eq. 7.7\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P\left( {{\mathrm{diff}}} \right) = \frac{{{\mathrm{tC}}^{\mathrm{n}}}}{{{\mathrm{tC}}^{\mathrm{n}} + {{K}}_{{\mathrm{th}}}^{\mathrm{n}}}}$$\end{document}Pdiff=tCntCn+Kthn

tC: Threshold count, number of timesteps a cell is over Cth\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{l}ie.{\mathrm{IF}}\,C > C_{\mathrm{th}}:{\mathrm{tC}} = {\mathrm{tC}} + 1\,{\mathrm{ELSE}}:{\mathrm{tC}} = {\mathrm{tC}} - 1\\ {{K}}_{{\mathrm{th}}} = 2\\ {{n}} = 4\end{array}$$\end{document}ie.IFC>Cth:tC=tC+1ELSE:tC=tC-1Kth=2n=4

A stochastic differentiation term is incorporated to represent external factors affecting differentiation that are not mechanistically captured by our intercellular diffusion module. Since the main external factor was the addition of retinoic acid, which was replaced every 24 h, the stochastic term is made oscillatory to correspond with RA-addition. Specifically, differentiation potential increases for 12 h after RA-addition then decreases back to the base differentiation potential, with this repeating every 24 h (Eq. 8).8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P\left( {{\mathrm{stochastic}}} \right) = 0.002 + {\mathrm{Abs}}\left[ {0.001 \ast {\mathrm{sin}}\left( {\frac{{t_{{\mathrm{sim}}}}}{8}} \right)} \right]$$\end{document}Pstochastic=0.002+Abs0.001*sintsim8where tsim is the simulation time.

To convert our experimental images into digital networks, we used a CellProfiler (http://cellprofiler.org/) pipeline (Supplementary Figures 1, 6). The pipeline uses the blue channel (Hoechst) to identify nuclei within the image, and then uses propagation to define cellular boundaries. Both red and green channels, representing Sox2 and Oct4, respectively, first use a local MCT threshold to separate high expressers from low expressers and background noise. After thresholding, the red and green channels are converted to masks and applied to the cell objects identified using the blue channel. This data is exported to a CVS file and imported into python, where it is represented in an analogous format as our simulation data.

We computationally generated 2D pattern classes defined for embryoid body differentiation. Each pattern class was produced using a previously digitized experimental network, resulting in 120 unique colony structures. To better distinguish spatially localized clusters, the differentiation status of a cell was also a function of the differentiation status of neighboring cells. Specifically, if the percent of differentiated neighbors is greater than the average number of differentiated neighbors for every cell in the population, it is classified as differentiated. From each network, seven metrics were extracted as represented in Supplementary Figure 6. After collecting these metrics for every pattern class generated in our training set (7 metrics × 960 networks), the data is mean centered, scaled to unit variance, then transformed using Principal Component Analysis (PCA). The same metrics were extracted from both simulation and experimental data, scaled/normalized in an identical fashion to our training set, and transformed into latent variable space using the previously trained PCA transform. The average trajectories for simulations were calculated by averaging every data point at a specific time.

All tests of significance, unless otherwise stated, were conducted via t-test with a significance level (α) of 0.01. Samples for the analysis of the differentiation trajectory and perturbation studies were collected from three different culture plates for each time and perturbation, with each plate contributing 10–20 images to the final n-values listed in their respective figure legends.

The computational model can be accessed from Github (DOI: 10.5281/zenodo.1413539) Software utilized: Python 2.7, CellProfiler 2.1.1

Supplementary Information Description of Additional Supplementary Files Supplementary Movie 1 Supplementary Movie 2