(C) PLOS One This story was originally published by PLOS One and is unaltered. . . . . . . . . . . A computational modeling approach for predicting multicell spheroid patterns based on signaling-induced differential adhesion [1] ['Nikita Sivakumar', 'Department Of Biomedical Engineering', 'University Of Virginia', 'Charlottesville', 'Virginia', 'United States Of America', 'Helen V. Warner', 'Shayn M. Peirce', 'Matthew J. Lazzara', 'Department Of Chemical Engineering'] Date: 2022-12 Physiological and pathological processes including embryogenesis and tumorigenesis rely on the ability of individual cells to work collectively to form multicell patterns. In these heterogeneous multicell systems, cell-cell signaling induces differential adhesion between cells that leads to tissue-level patterning. However, the sensitivity of pattern formation to changes in the strengths of signaling or cell adhesion processes is not well understood. Prior work has explored these issues using synthetically engineered heterogeneous multicell spheroid systems, in which cell subpopulations engage in bidirectional intercellular signaling to regulate the expression of different cadherins. While engineered cell systems provide excellent experimental tools to observe pattern formation in cell populations, computational models of these systems may be leveraged to explore more systematically how specific combinations of signaling and adhesion parameters can drive the emergence of unique patterns. We developed and validated two- and three-dimensional agent-based models (ABMs) of spheroid patterning for previously described cells engineered with a bidirectional signaling circuit that regulates N- and P-cadherin expression. Systematic exploration of model predictions, some of which were experimentally validated, revealed how cell seeding parameters, the order of signaling events, probabilities of induced cadherin expression, and homotypic adhesion strengths affect pattern formation. Unsupervised clustering was also used to map combinations of signaling and adhesion parameters to these unique spheroid patterns predicted by the ABM. Finally, we demonstrated how the model may be deployed to design new synthetic cell signaling circuits based on a desired final multicell pattern. The remarkable ability of cells to self-organize is critical for the assembly of functional tissues during embryogenesis and is impaired by the molecular aberrations that lead to tumorigenesis and metastasis. Learning the rules of cellular self-assembly will provide a new way to understand such physiological and pathophysiological processes and create an instruction manual for designing tissues from scratch for therapeutic applications. Biologists have begun to learn the rules of self-assembly through reverse engineering–that is, through engineering biochemical circuits that control how cells adhere to one another. These approaches can yield simple, multi-cell structures that self-assemble into a core of one cell type surrounded by a shell of another cell type, for example. However, engineering more complex tissue patterns requires exploring a large domain of circuit structures and parameters. To facilitate this exploration in a guided and systematic manner, we created a computational model that predicts how mixtures of cells with different circuity for expressing adhesion proteins will interact to form varied patterns in spheroidal tissues comprised of hundreds of cells. Our model enables the design of internal molecular signaling circuitry that permits cells to be used as building blocks for self-assembled tissues with specific structures, and ultimately functions. Funding: This work was supported by funding from the National Science Foundation Grant No. 1700687 (MJL), University of Virginia Engineering in Medicine Program (MJL, SPC), University of Virginia Center for Advanced Biomanufacturing (MJL), and Arnold and Mabel Beckman Foundation (NS). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Copyright: © 2022 Sivakumar et al. This is an open access article distributed under the terms of the Creative Commons Attribution License , which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. While it would be impossible to explore the effects of perturbing different system parameters exhaustively through experiments, computational models provide an efficient method for predicting unique multicellular patterns for different parameter combinations [ 42 – 46 ]. ABMs are a particularly useful computational approach for modeling heterogeneous cell types and spheroid patterning because ABMs can represent individual cells within a population and predict how cell type-specific behaviors produce spatial patterns at the multicell level [ 8 – 10 , 28 , 47 – 54 ]. Here, we developed two-dimensional (2D) and three-dimensional (3D) ABMs of the Ncad/Pcad bidirectional signaling system engineered by Toda et al. [ 2 ]. Cells were modeled as agents expressing different cadherins with a specified time delay in response to heterotypic receptor-ligand interactions and as differentially adhering to like or unlike cells based on the homotypic and heterotypic adhesion properties of the cadherins. Model parameters were tuned to fit the ABM to experimental observations we made for a 1:1 ratio of the seeded cell types, and the fitted model was experimentally validated against experiments with unequal cell seeding densities. Unsupervised machine learning was then employed to map different combinations of perturbed model parameters (e.g., characterizing cadherin induction and adhesion strengths) to specific spheroid patterns predicted by the ABM, including the previously described core/pole and core/shell patterns and others resembling a soccer ball, bull’s eye, or stripes. The ABM was also used to design synthetic cell-cell signaling circuits that encode specified multicellular patterns, demonstrating a potential practical application of this modeling approach. (A) A schematic of the synthetic bidirectional signaling circuit engineered by Toda et al. [ 2 ] and implemented in the agent-based model (ABM) is shown. Spheroids are seeded with non-adhesive “sender” BFP+ cells and “receiver” -/- cells (labeled -/- to indicate no fluorescent protein expression). Sender BFP+ cells can induce GFP/Ncad expression in -/- receiver cells, and newly formed GFP/Ncad cells behave as sender cells to induce mCherry/Pcad expression reciprocally in BFP+ cells (at which point BFP+ cells behave as receiver cells). The final spheroid composition can include GFP/Ncad, mCherry/Pcad, or BFP+ cells (that never contacted a GFP/Ncad cell and thus failed to transition to the mCherry/Pcad phenotype). In this ruleset, the model cadherin induction constant (c i ) controls the probability of GFP/Ncad induction in -/- cells based on neighboring BFP+ cells and the probability of mCherry/Pcad induction in BFP+ cells based on neighboring GFP/Ncad cells. The model homotypic adhesion strength (P i ) controls the probability of forming junctions between cells expressing like cadherins on adjacent grid spaces. (B) Metrics used to analyze ABM outputs are represented schematically. Regions of a given color comprising three or more cells are defined as cell clusters. These metrics are further clarified in Table 1 . (C) The rules executed at each ABM timestep are depicted. Parameters, indicated in bold italics, are defined in Table 2 . Synthetic biology approaches have been used to impart self-organizing properties to multicellular systems [ 2 , 35 – 41 ]. These systems provide an ideal experimental setting to explore pattern formation resulting from dynamic cell-cell signaling because the cell interactions are clearly defined and programmable. A recent example of interest is the bidirectional signaling system engineered by Toda et al. [ 2 ], in which sender cells drive expression of a particular cadherin in receiver cells that then reciprocally drive expression of a different cadherin in the initial sender cells ( Fig 1A ). The authors observed that the relative homotypic and heterotypic cadherin interaction strengths directly impacted observed spheroid patterns. When the network was engineered with N-cadherin (Ncad) and P-cadherin (Pcad), which both have high homotypic and low heterotypic adhesion strengths, a “core/pole” pattern arose. Specifically, Ncad-expressing receiver cells formed a core decorated by smaller clusters, or poles, of Pcad-expressing sender cells and surrounding non-adhesive sender cells [ 2 ]. When the system was instead engineered with high and low expression of E-cadherin (Ecad), creating cells with high (but differential) homotypic and heterotypic adhesion, a “core/shell” pattern arose. Specifically, high-Ecad receiver cells formed a core encapsulated by a shell of low-Ecad sender cells [ 2 ]. In both systems, the number of poles and core/shell spheroids formed varied across trials with the same initial conditions [ 2 ], but the factors that led to this variable outcome were not systematically investigated. These results open several questions about how combinations of cell-cell signaling and differential adhesion events affect multicell patterning. The multicell patterns that form in these systems may depend on the initial ratio of cell subpopulations, the order of signaling events, and/or the induced adhesion strengths between individual cells, but more work is needed to explore such potential dependencies and how combinations of these perturbations may affect the robustness of pattern formation. The ability of a cell to adhere to another cell may be time- and space-dependent based on cell-cell signaling that alters expression of adhesion proteins. In embryonic development, for example, cell-cell interactions initiate transcriptional processes that regulate expression of specific adhesion proteins, leading to the segregation of cell subpopulations into organ-specific tissue layers [ 12 – 19 ]. Similarly, in many cancers, signaling cascades cause individual cells to downregulate expression of the cell adhesion protein E-cadherin, leading to epithelial-mesenchymal transition (EMT) [ 20 – 24 ]. EMT creates the opportunity for differential adhesion between epithelial cells and the mesenchymal counterparts they can generate, impacting the spatial configuration of these subpopulations throughout the tumor [ 25 – 28 ]. Emergent cell patterning within tumors can impact tumor growth, migration, and response to treatment [ 20 , 28 – 34 ]. Understanding the mechanisms by which multicell patterning occurs can inform drug targets for various congenital diseases and other pathologies that occur throughout adulthood. Multicell patterning processes in embryonic development and tumorigenesis are regulated by dynamic cell-cell signaling and differential adhesion between cell subpopulations [ 1 , 2 ]. Steinberg and colleagues [ 3 ] previously proposed the differential adhesion hypothesis to explain how patterns emerge at the multicell level in heterogenous mixtures of cells. Specifically, they proposed that cell subpopulations self-organize to minimize the adhesive free energy between cells such that more adhesive cell subpopulations are surrounded by less adhesive ones [ 3 – 5 ]. The differential adhesion hypothesis has been validated experimentally [ 3 – 7 ] and incorporated into computational models of cell sorting and morphogenesis [ 8 – 11 ]. However, this conceptual model does not account for how dynamic cell-cell signaling processes, which regulate adhesion protein expression, may yield emergent patterns that differ from those predicted by purely equilibrium considerations (i.e., minimization of adhesive free energy). Results 2D and 3D ABMs predict the impacts of dynamically activated differential adhesion on spheroid patterning 2D and 3D ABMs were developed based on the synthetic Pcad/Ncad circuit engineered by Toda et al. [2] (Fig 1A). The models were seeded with “sender” and “receiver” cells. Sender cells express blue fluorescent protein (BFP+), CD19 ligand, and green fluorescent protein (GFP) receptor. Receiver cells express the CD19 receptor (-/-). In response to CD19 receptor-ligand binding between the two cell types, receiver -/- cells express Ncad and GFP to become GFP/Ncad cells. GFP/Ncad cells then behave as sender cells and reciprocally promote Pcad and mCherry expression in BFP+ cells (that then behave as receiver cells) via GFP ligand-receptor binding, forming mCherry/Pcad cells. Hence, in this synthetically engineered multicell system, both cell types can function as either sender or receiver cells based on receptor-ligand binding interactions with adjacent cells, and the terms “sender” and “receiver” are used to clarify a cell’s role during specific steps in bidirectional signaling. Final spheroid patterns were quantified using metrics such as total fraction of spheroid area occupied by each cell type, average areas of homotypic cell clusters, and average distances between homotypic cell clusters and the spheroid center (Fig 1B, Table 1). Homotypic cell “clusters” were defined as three or more cells of the same cell type adjoined by cell-cell junctions on adjacent grid spaces in the cardinal directions. The “core” cell type of a spheroid was defined as the cell type of the homotypic cluster with minimum distance to all other homotypic clusters in the spheroid. PPT PowerPoint slide PNG larger image TIFF original image Download: Table 1. Metrics used to quantify multicell patterns in 2D ABMs. https://doi.org/10.1371/journal.pcbi.1010701.t001 Interactions among cells in the ABM were encoded as six rules incorporating nine model parameters that execute on each timestep of the simulation, where one timestep represents 30 min (Fig 1C, Table 2). Cadherin phenotypes were induced in the model based on induction constants (c i ) that control the probability of cadherin induction in a receiver cell based on the number of neighboring sender cells. Over the time course of the model, neighboring cadherin-expressing cells form junctions with each other based on homotypic (P i ) and heterotypic (P a,b ) adhesion strengths encoded as probabilities of adhesion in the model. To recapitulate the Ncad/Pcad spheroids observed in Toda et al. [2], we assumed that both cadherin cell types have an equal induction constant (c i ) and equal homotypic adhesion strength (P i ). In later sections, we explore how spheroid patterning is impacted by assigning unequal values of c i and P i to the GFP/Ncad phenotype (c a , P a ) and to the mCherry/Pcad phenotype (c b , P b ). All simulations that were compared to in vitro experiments were run for 100 timesteps, where one timestep represents 30 min, because in vitro spheroids were observed for 50 hr. Additionally, because published data in Toda et al. [2] indicated that cells roughly doubled in number over the course of 50 hr, ABMs were seeded with twice the number of cell agents as the corresponding in vitro experiment (i.e., an in vitro experiment seeded with 200 cells corresponded to an ABM seeded with 400 cell agents). All other simulations were run based on pilot simulations to identify runtimes beyond which no changes were observed in the spheroid pattern, and the specific model runtime for each set of simulations is denoted in corresponding figure captions. PPT PowerPoint slide PNG larger image TIFF original image Download: Table 2. ABM parameter values. https://doi.org/10.1371/journal.pcbi.1010701.t002 Published data and an evolutionary algorithm were used to parameterize the ABM Parameters in the 2D ABM were either estimated based on published experimental data in Toda et al. [2] or computationally tuned with an evolutionary algorithm (EA) (Table 2). The time delay for cadherin phenotypes to be induced (a for GFP/Ncad and b for mCherry/Pcad) and the times to express each cadherin (Δa, Δb) were manually adjusted in the 2D ABM according to published data in Toda et al. [2]. Specifically, these parameters were set such that in simulations of experiments seeded with 100 BFP+ and 100 -/- cells, GFP/Ncad cells formed after 13 hr and mCherry/Pcad cells developed adjacent to GFP/Ncad cells after 21 hr, the respective times reported by Toda et al. [2]. Because the remaining five parameters (c i , P i , P a,b , μ, t) could not be determined from experimental data, these parameters were first manually approximated in the 2D ABM as values yielding outputs qualitatively similar to published images of spheroids from the Ncad/Pcad bidirectional signaling system in Toda et al. [2] (i.e., characterized by a GFP/Ncad central core flanked by mCherry/Pcad clusters and individual BFP+ cells). Then, a univariate parameter sensitivity analysis was implemented for the 2D ABM to determine which of these parameters had the largest relative impact on spheroid patterning. Each of the five parameters was increased or decreased by 10% from its base value, and the effects of these perturbations on the numbers and average areas of green and red clusters were determined (Fig 2A–2D). This analysis identified the cadherin induction constant (c i ) and the homotypic adhesion strength of both cadherin cell types (P i ) as the two parameters whose variation produced the largest changes in all four measured outputs. In this univariate sensitivity analysis, we assumed both cadherin cell types have an equal induction constant (c i ) and equal homotypic adhesion strength (P i ). PPT PowerPoint slide PNG larger image TIFF original image Download: Fig 2. Parameter sensitivity analysis and optimization improves model fit to in vitro data. Parameters that could not be estimated based on experimental data (c i , P i , P a,b , μ, t) were perturbed by a factor 0.9 and 1.1, and the sensitivity coefficient was quantified for the following outputs: (A) number of green clusters, (B) average fractional area of green clusters, (C) number of red clusters, and (D) average fractional area of red clusters. Sensitivity coefficients are reported as the average of 100 model runs, each run for 100 timesteps. The x-axis is ordered from greatest to least absolute-value sensitivity coefficient for each panel. (E) An evolutionary algorithm (EA) was used to tune parameters c i and P i by minimizing the error between model outputs and experimental images of 200-cell spheroids seeded with a 1:1 ratio of BFP+ and -/- cells. The EA was run 30 times for a population of 20 model combinations over 100 generations. Within each generation, the error for each parameter combination was averaged over 10 runs for 100 timesteps. Sample model graphical outputs generated from the parameter combinations with maximum or minimum error are shown. https://doi.org/10.1371/journal.pcbi.1010701.g002 Values for c i and P i were determined using an evolutionary algorithm (EA) that optimized the 2D ABM prediction of patterns that were experimentally observed in spheroids seeded with 200 cells at a 1:1 ratio of BFP+ and -/- cells. A detailed description of the EA implementation is provided in Materials and Methods. A plot of 30 EA runs over 100 generations demonstrates that the EA accomplished a substantial reduction in error (difference between model prediction and data) to identify values for c i and P i (Fig 2E). The parameter values optimized for the 2D ABM were also used in the 3D ABM. Dynamic cell-cell signaling and homotypic adhesion strength impact pattern formation To test the hypothesis that the dynamic expression of different cadherin isoforms impacts spheroid patterning, three different signaling circuits were explored (Fig 6A). In circuit 1, the cadherins were constitutively expressed, as opposed to expression via induction by juxtracrine signaling. In circuit 2, only unidirectional signaling was present such that seeded mCherry/Pcad sender cells induced the GFP/Ncad phenotype in seeded -/- receiver cells. In circuit 3, the full bidirectional signaling circuit was modeled, with seeded BFP+ sender cells inducing GFP/Ncad expression in seeded -/- receiver cells and newly formed GFP/Ncad cells inducing mCherry/Pcad expression in BFP+ cells. In this circuit, the induction constants for GFP/Ncad and mCherry/Pcad phenotypes were parametrized separately (c a and c b , respectively, instead of c i ). In the simulations discussed for circuit 3 in this section, c a = c b = c i , but c a and c b are varied in later Results sections. For each circuit, three homotypic adhesion strength profiles were tested, in which the homotypic adhesion strengths of GFP/Ncad and mCherry/Pcad phenotypes were parametrized separately (P a and P b , respectively, instead of P i ). In profile A, GFP/Ncad and mCherry/Pcad cell types exhibited maximal homotypic adhesion and minimal heterotypic adhesion (P a = P b = 1.0, P a,b = 0.0), giving rise to core/pole spheroids, in which one adhesive cell type assembled in a core and the other assembled in smaller peripheral clusters. In profile B, GFP/Ncad cells exhibited maximal homotypic adhesion, mCherry/Pcad cells exhibited minimal homotypic adhesion, and heterotypic adhesion was maximal (P a = 1.0, P b = 0.0, P a,b = 1.0). In profile C, GFP/Ncad cells exhibited minimal homotypic adhesion, mCherry/Pcad cells exhibited maximal homotypic adhesion, and heterotypic adhesion was maximal (P a = 0.0, P b = 1.0, P a,b = 1.0). Both profiles B and C gave rise to core/shell spheroids, in which the more adhesive cell type formed a core surrounded by a complete shell of the less adhesive cells (Fig 6A). Combinations of each signaling circuit and adhesion profile yielded 9 rulesets. For all rulesets, simulations were seeded with a 1:1 ratio of initial cell types for a fixed total of 200 cells. PPT PowerPoint slide PNG larger image TIFF original image Download: Fig 6. Dynamic cell-cell signaling has differential effects on core/pole and core/shell spheroids. (A) A schematic of three different ABM signaling circuits and three adhesion profiles is shown. Circuit 1 contains constitutive cadherin expression and no juxtracrine signaling. Circuit 2 includes unidirectional signaling. Circuit 3 includes bidirectional signaling (as described in Fig 1A). In Circuit 3, induction constants for GFP/Ncad and mCherry/Pcad phenotypes are parametrized separately (c a and c b , respectively). In each adhesion profile, the homotypic adhesion strengths of GFP/Ncad and mCherry/Pcad phenotypes are also parametrized separately (P a and P b , respectively). All depicted adhesion profiles have a homotypic adhesion strength of 1.0. Adhesion profile A gives rise to core/pole spheroids, whereas adhesion profiles B and C give rise to core/shell spheroids. The combination of circuit type and adhesion profile leads to nine rulesets (1A, 1B, 1C, 2A, 2B, 2C, 3A, 3B, and 3C). (B) The model was run 100 times, each for 100 timesteps, for each of the nine rulesets. Each ruleset was seeded with 200 cells at a 1:1 ratio of the seeding cell types. The fraction of trials with either GFP/Ncad or mCherry/Pcad spheroid core was compared between rulesets using a Chi-Squared test with df = 8; significant differences between observed and expected fraction (0.50) of core cell types were found across all nine rulesets with p < 0.001. The expected fraction was set to 0.5 because the null hypothesis states that signaling has no effect on which cell type forms the spheroid core. * indicates groups with residuals > 1. https://doi.org/10.1371/journal.pcbi.1010701.g006 The effects of uni- and bi-directional signaling between cells were explored by comparing the core cell types in spheroids for circuits 1–3 (Fig 6B). Adhesion profile A tended to produce core/pole spheroids for which the core cell type varied based on the signaling circuit. Ruleset 1A gave rise to either GFP/Ncad or mCherry/Pcad cores, whereas in Rulesets 2A and 3A the adhesive phenotype that was first present in the circuit formed the spheroid core. Thus, the adhesive cell type that is present longer in the spheroid exhibited a higher likelihood of condensing into a central core. The same trends for core cell type were not observed for adhesion profiles B and C, which tended to produce core/shell patterns. In core/shell spheroids, the more adhesive phenotype formed the spheroid core in > 90% of model simulations regardless of the signaling circuit (Fig 6B). Thus, for profiles B and C, the signaling dynamics produced outcomes that were consistent with the differential adhesion hypothesis in 200-cell spheroids seeded with a 1:1 ratio. Varying cell seeding in Ruleset 3B leads to patterns that violate the differential adhesion hypothesis To test whether it is possible to generate spheroids that violate the differential adhesion hypothesis, the 2D ABM was used to predict patterns of cell subpopulations when the seeding ratios and total cell numbers in Ruleset 3B were varied. Ruleset 3B was chosen because it is based on the fundamental premise of the differential adhesion hypothesis, in which two phenotypes with differing homotypic adhesion strengths self-organize. As the relative proportion of seeded -/- cells and spheroid radius increased, the likelihood that adhesive cells would enclose non-adhesive cells increased (Fig 7A). In the 1:1 ratio group, regardless of spheroid radius, a small number of simulated spheroids (~12–20%) contained one group of non-adhesive red cells enveloped by adhesive green cells, while virtually no spheroids contained non-adhesive -/-, or gray, cells enveloped by adhesive green cells. For the 1:4 and 1:9 ratios, the fraction of spheroids with groups of non-adhesive red cells enclosed by adhesive green cells decreased significantly, and the overall fraction of spheroids with groups of non-adhesive gray cells enclosed by adhesive green cells increased significantly, with the latter effect becoming more apparent as the spheroid radius increased. This trend is partially explained by the presence of fewer seeded BFP+ cells, which leads to fewer -/- cells transitioning to the GFP/Ncad phenotype and mCherry/Pcad phenotypes. Due to the increased tendency of smaller green clusters to spontaneously arise at the periphery of the spheroid, non-adhesive cells surrounding each individual green cluster appeared to be trapped in the center of the spheroid, giving rise to situations where non-adhesive cells appeared enclosed by adhesive cells. This interpretation is reinforced by the fact that when the induction constant of the GFP/Ncad phenotype (c a ) was increased, a single, larger green core formed, and the spheroid patterns resembled those observed in the 1:1 group (Fig 7B and 7C). When the induction constant of the GFP/Ncad phenotype (c a ) was larger, fewer BFP+ sender cells were required to induce the GFP/Ncad phenotype in -/- receiver cells, and a larger central GFP/Ncad core formed more quickly in the spheroid, instead of the slower formation of GFP/Ncad clusters and their subsequent enclosure of non-adhesive cells. PPT PowerPoint slide PNG larger image TIFF original image Download: Fig 7. Increasing the relative fraction of seeded -/- cells in Ruleset 3B leads to instances where non-adhesive cell types are surrounded by adhesive cell types. (A) The number of gray and red non-adhesive cell clusters with > 75% enclosure by green adhesive cells was averaged for 100 model runs, each run for 200 timesteps, for spheroids with radii 5.6, 6.4, and 7.2 pixels (px) and varying initial BFP+:-/- ratios (1:9, 1:4, 1:1, 4:1, 9:1). All groups were compared with a Chi-Squared test with df = 14, resulting in a p < 0.001; * indicates groups with greater than 10% contribution to the chi statistic. (B) c a and c b were simultaneously varied over values of 0.05, 0.2, and 10 in spheroids with radius 7.2 px and initial BFP+:-/- ratio of 1:4. Total number of gray clusters with > 75% enclosure by green cells was quantified for 100 model runs, each run for 200 timesteps, for each model condition. (C) A representative output of a simulated spheroid with radius 7.2 px, initial BFP+:-/- seeding ratio of 1:4, and c a = 10 and c b = 0.05. https://doi.org/10.1371/journal.pcbi.1010701.g007 Homotypic adhesion strengths impact heterogeneity in spheroid patterns Systematically covarying the homotypic adhesion strengths of GFP/Ncad and mCherry/Pcad cells (P a and P b , respectively) and setting P a,b = 0 in the 2D ABM for signaling circuit 3 (bidirectional signaling) for spheroids seeded with ~350 cells produced a range of different spheroid patterns. We chose to increase the spheroid size for these models because distinct patterns emerged more clearly in larger spheroids. To display the results of these simulations, heatmaps were created to represent the numbers and average areas of green or red clusters in the spheroids (Fig 8A). Recall that the extremes of these maps describe core/shell (P a = 1.0 and P b = 0.0 or P a = 0.0 and P b = 1.0) and core/pole patterns (P a = P b = 1.0). With P a > P b , there were typically fewer and larger green clusters than red clusters, but the average number and area of green clusters were not as low as those observed in the core/pole pattern with a red core. Conversely, with P b > P a , there were fewer and larger red clusters than green clusters, but the average number and area of red clusters were not as low as those observed in the core/pole pattern with a green core. Deviations from core/pole and core/shell patterns could occur when P a and P b have non-extreme values. PPT PowerPoint slide PNG larger image TIFF original image Download: Fig 8. Covarying homotypic adhesion strengths of activated cell types in Ruleset 3A leads to a variety of different patterns. (A) Homotypic adhesion strengths of the first and second activated phenotypes (P a , P b ) were simultaneously varied with values of 0.0, 0.25, 0.50, 0.75, and 1.0 for signaling circuit 3 (bidirectional signaling circuit) in the 2D ABM. For each parameter combination, the simulations were run for spheroid radius 6.4 px at a 1:1 BFP+:-/- ratio and were run 100 times, each for 300 timesteps. The average values of the indicated metrics were plotted in heat maps. (B) Results from this parameter space exploration were clustered using a k-means algorithm into six groups, with representative ABM outputs shown for each group and color coded as the key for this graph. (C) A t-distributed stochastic neighbor embedding (tSNE) algorithm was used to show where the k-means clusters projected in two dimensions. The tSNE algorithm was run with the perplexity parameter = 100. https://doi.org/10.1371/journal.pcbi.1010701.g008 To identify the different spheroid patterns formed, we performed k-means clustering, an unsupervised machine learning technique, to identify six groups within the data (Fig 8B). Each of the six groups corresponded to a distinct spheroid pattern: a core/pole pattern with a green core and three red poles, a core/shell pattern with green core and red shell, a core/shell pattern with red core and green shell, a striped pattern of alternating green and red linear clusters, a soccer ball pattern of circular green and red clusters, and a bull’s eye pattern with radially alternating green and red shells. As an alternative way to visualize the k-means clustering results, the data were projected in two dimensions using t-distributed stochastic neighbor embedding (tSNE), with pattern identities retained from k-means clustering (Fig 8C). These results also visually support the binning of spheroids into more than just two categories. Each combination of homotypic adhesion strengths led to simulation outputs that exhibited varying frequencies of the six patterns identified by the clustering algorithm described above. When P a = P b = 0.0, a soccer ball pattern was generated more than 90% of the time. In alignment with the differential adhesion hypothesis, when P b > P a and P a = 0.0, patterns with a red core were generated almost 100% of the time. Conversely, when P a > P b and P a = 0.0, patterns with a green core were generated almost 100% of the time. The distribution of core/shell, core/pole, bull’s eye, and striped patterns generated was relatively similar for all conditions where the homotypic adhesion strengths of both cell types are equal (e.g., [P a , P b ] = [0.25, 0.25], [0.5, 0.5], [0.75, 0.75], and [1.0, 1.0]). Interestingly, when P a = P b = 1.0, patterns with a green core and surrounding red cells formed, but a bull’s eye pattern emerged with equal frequency. In larger spheroids, the possibility for there to be alternating clusters of cell types was higher. When the GFP/Ncad phenotype had relatively weak homotypic adhesion (e.g., [P a , P b ] = [0.25, 0.25] or [0.25, 0.75]), striped and bull’s eye patterns emerged more frequently than in conditions when the GFP/Ncad phenotype had stronger homotypic adhesion (e.g., [P a , P b ] = [0.75, 0.25] and [0.75, 0.75]). When the GFP/Ncad phenotype had relatively weak homotypic adhesion, the frequency of heterotypic interactions between cells was greater, allowing for a higher frequency of patterns characterized by longer borders between heterotypic cells. In conditions where differential adhesion was present (e.g., [P a , P b ] = [0.25, 0.75] and [0.75, 0.25]), the most frequent patterns involved the more adhesive cells being surrounded by the less adhesive cells. With P a = 0.75 and P b = 0.25, the most frequent patterns were the core/pole and core/shell spheroids with a green core and bull’s eye spheroids with a green core. With P b = 0.75 and P a = 0.25, the most frequent patterns were the core/pole spheroids with a red core and the bull’s eye spheroids with a red core. This trend suggests that when differential adhesion is present, but the adhesive strength of the second induced phenotype is not maximal, hybrids of the core/shell pattern arise. Cadherin induction constants impact core/pole and core/shell patterns Varying the induction constants of GFP/Ncad and mCherry/Pcad cells (c a and c b , respectively) was also predicted to impact the final numbers of green and red clusters and their average areas (Fig 9A). For Ruleset 3A (bidirectional signaling core/pole ruleset) and when c a = 0.05 and c b = 10, one green cluster with a relatively small area was formed, and one red cluster with a larger area was formed. When c b was only slightly larger than c a , multiple small green clusters were formed, and one or two red clusters were formed, suggesting the possibility that different patterns may occur for different parameter combinations. Seven groups were identified by k-means clustering (Fig 9B). The first two patterns included a green core and three or four red poles. The third pattern comprised spheroids with a green core and red shell. The fourth pattern exhibited an inversion of the first two, with a red core and two or more green poles. The fifth pattern included one red core and one green pole. The sixth included a red core and single green and gray peripheral cells. Finally, the seventh pattern included striped or soccer ball patterns. As in Fig 9A, when c a << c b , a red core formed with one or multiple green poles. However, when c a was sufficiently high, core/pole structures with a green core and multiple red poles formed more frequently. To produce an alternative visualization, tSNE was used to reduce the dimensionality of the system (Fig 9C). tSNE results further supported the existence of multiple subtypes of spheroids. Modifying the cadherin induction constants affects the overall spheroid pattern because various combinations of c a and c b directly impact the final ratio of BFP+, -/-, GFP/Ncad, and mCherry/Pcad cells in the spheroid (S1A Fig). PPT PowerPoint slide PNG larger image TIFF original image Download: Fig 9. Cadherin induction constants (c i ) control the number of poles in core/pole spheroids. (A) The cadherin induction constants (c a , c b ) were simultaneously varied across values of 0.05, 0.2, and 10 for Ruleset 3A (bidirectional signaling core/pole ruleset) in the 2D ABM. For each parameter combination, simulations were seeded with 200 cells at a 1:1 BFP+:-/- ratio and run 100 times, each for 100 timesteps. The average from 100 model runs of the following metrics were plotted for each parameter combination: number of green clusters, number of red clusters, average fractional area of green clusters, and average area of red clusters. (B) Results from this parameter space exploration were clustered using a k-means algorithm into seven groups, with representative ABM outputs shown for each group and color coded as the key for this graph. (C) t-distributed stochastic neighbor embedding (tSNE) algorithm was used to project the data into two dimensions with colors for spheroid cluster types retained from the k-means clustering analysis in panel (B). The tSNE algorithm was run with perplexity = 100. https://doi.org/10.1371/journal.pcbi.1010701.g009 The parameters c a and c b were also varied for Ruleset 3B (bidirectional signaling core/shell ruleset). The ratio of the shell’s thickness to the core’s radius progressively increased as c b became larger than c a (S1B and S1C Fig), which is consistent with more mCherry/Pcad cells at the final time point. Because this parameter describes sensitivity of a cell to its neighbors, and the number of neighbors in 3D is larger than in 2D, results from these 2D model calculations may not translate directly to a 3D model. [END] --- [1] Url: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1010701 Published and (C) by PLOS One Content appears here under this condition or license: Creative Commons - Attribution BY 4.0. via Magical.Fish Gopher News Feeds: gopher://magical.fish/1/feeds/news/plosone/