Continuous dynamic modeling of regulated cell adhesion

Cell-cell adhesion is essential for tissue growth and multicellular pattern formation, and crucial for the cellular dynamics during embryogenesis and cancer progression. Understanding the dynamical gene regulation of cell adhesion molecules (CAMs) responsible for the emerging spatial tissue behaviors is a current challenge due to the complexity of these non-linear interactions and feedback loops at different levels of abstraction—from genetic regulation to whole-organism shape formation. Continuous mathematical models of cell adhesion are ideal for the modeling of the spatial dynamics of large cell populations, where different cell types define inherent adhesion strengths. However, biologically the adhesive properties of the cell arise dynamically from differential expression of CAMs, which are precisely regulated during development and cancer progression. To extend our understanding of cell and tissue behaviors due to the regulation of adhesion molecules, here we present a novel model for the spatial dynamics of cellular patterning, growth, and shape formation due to the differential expression of CAMs and their regulation. Capturing the dynamic interplay between genetic regulation, CAM expression, and differential cell adhesion, the proposed continuous model can recapitulate the complex and emergent spatial behaviors of cell populations that change their adhesion properties dynamically due to inter- and intracellular genetic regulation. This approach can demonstrate the mechanisms responsible for classical cell sorting behaviors, cell intercalation in proliferating populations, and the involution of germ layer cells induced by a diffusing morphogen during gastrulation. Integrating the emergent spatial tissue behaviors with the regulation of genes responsible for essential cellular properties such as adhesion will pave the way towards understanding the genetic regulation of large-scale complex patterns and shapes formation in developmental, regenerative, and cancer biology.


Introduction
The adhesive properties of cells can dictate their spatial behaviors and the formation of correct tissue patterns and shapes during morphogenesis and homeostasis (1). Seminal studies demonstrated how stirred disassociated embryonic tissues could sort themselves and regain their specific configurations due to the distinct adhesive properties of their different cell types (2,3).
These cell-cell adhesive forces are dependent on the expression of cell adhesion molecules (CAMs) through the cell surface, such as families of proteins including the cadherins, integrins, and nectins (4,5). CAMs expressed at the cell surface can form bonds with the same or different CAM types expressed in neighbor cells, resulting in different adhesive strengths. These CAM adhesive forces are transmitted to the cell through its cytoskeleton network and can result in specific cell spatial behaviors. The sum of intercellular interactions between different CAMs determine the net force in the cell, which drive specific cellular movements and emergent tissue patterns. The importance of cell adhesion is clear when its cellular components are perturbed, resulting in tissues that can degenerate into mis-patterned phenotypes during development (6) and disease states such as cancer progression and metastasis (7,8). However, the biophysical dynamics and cellular behaviors directed by differential adhesion and its genetic regulation are currently not completely understood.
The precise regulation of CAM expression modulates the adhesive properties of cells and hence can control the movement of cells and the formation of global tissue patterns during morphogenesis, whereas its dysregulation may lead to tumor formation and metastasis. Several gene families have been found to regulate CAM expression. The Snail family of transcription factors regulate the expression of cadherins essential for gastrulation in invertebrates, the epithelial-to-mesenchymal transition in neural crest cells in all amniotes, and the development of organs such as the kidney (9,10). Differential regulation of CAMs such as cadherins by ephrins and Hox genes is a key factor for proper cell distribution during limb morphogenesis and regeneration (11); mutations in these pathways can result in limbs with abnormal morphological organizations (12). Dysregulated pathways controlling CAMs expression is sufficient to induce tumor progression, metastasis formation, and drug resistance (9,13). Kinases can up-regulate Eselectin-a CAM essential for the localization of metastatic cancer cells in the lungs (14)-and specific kinase inhibitors targeting these pathways represent promising drugs for anticancer therapeutics (15). However, the complex feedback loops between CAM regulation, cellular adhesion dynamics, tissue behaviors, and intercellular signaling represents an extraordinary challenge that remains to be deciphered.
To understand the complex dynamics between the regulation of CAMs and the spatial tissue behaviors, mathematical and computational approaches are needed to model the physical properties of these processes and explain their emergent dynamics. Discrete models based on the extended Potts approach have been proposed to understand cell adhesion dynamics, and they can recapitulate the classical cell sorting dynamics due to adhesion (16)(17)(18), specific developmental dynamics (19)(20)(21), and cancer behaviors (22,23). These models do not include the dynamics of CAMs expression, and instead use pre-defined adhesion constants for different cell types.
Extensions to these discrete approaches have been proposed to model the concentration of CAMs, either using static concentrations defining cellular adhesion strengths (24) or dynamic concentrations with hybrid models (25). These approaches are based on the explicit modeling of cells, and hence computationally expensive for large numbers, which limits their applicability. In addition, mathematical methods to analyze discrete models are limited.
To overcome the limitations of discrete models, continuous models of cell adhesion have been proposed that can equally recapitulate the classical cell sorting behaviors but are computationally more efficient for the simulation of large populations and amenable for mathematical analysis (26,27). Continuous models have been successfully used to explain developmental processes (28) and cancer dynamics (29)(30)(31)(32). However, the adhesion properties in these models are static and defined with specific constants in pre-defined cell types. As a consequence, the regulation and dynamics of adhesion molecules have not been possible to model with continuous approaches, limiting our ability to understand the regulatory dynamics of CAMs expression and their influence in large scale tissue behaviors such as whole embryos.
Here we present a novel continuous model of cell adhesion due to the expression of CAMs and their regulation. This approach does not rely on pre-defined adhesion constants between cell types, but models as continuous the levels of CAM concentration, which in turn dynamically determine the adhesive properties between cells. Modeling the expression of CAMs naturally allows the inclusion of their regulatory dynamics, essential in many biological processes. We demonstrate the capabilities of the proposed model with three experiments. First, we show how the model can correctly recapitulate the classical Steinberg cell sorting dynamics due to the differential expression of CAMs. Next, we present a model of cellular intercalation dynamics resulting from the differential expression of two different nectins in a proliferating cell population. Finally, we model whole-embryo developmental dynamics during zebrafish gastrulation, showing how the diffusing morphogen Nodal regulates the expression of a cadherin, dynamically modulating the adhesive properties of cells and resulting in a characteristic involution of the mesendodermal germ layer. Integrating the regulatory dynamics of CAMs and their cell adhesion properties in the proposed continuous model permits the simulation and spatial predictions of the behaviors of large population of cells due to the interdependent dynamics of genetic regulation and adhesion proteins.

Model derivation
The dynamic adhesive properties of cells originate from the regulation and expression of CAMs.
CAMs expressed on neighbor cells interact with each other, generating adhesive forces. CAMs bind to both CAMs of the same type, as well as CAMs of other types. The adhesive force generated from interactions between CAMs hence depends on both the adhesive strength between CAMs and their specific levels of expression in the interacting cells. The dynamic regulation of CAM expression, possibly by intra-and intercellular regulatory factors, results in dynamic adhesive forces. These dynamic forces can dictate cellular and tissue movement, resulting in target patterns and shapes. The proposed model follows a continuous approach to define a population of cells with adhesion forces (26,27). The model does not explicitly define cell types, although it could, but instead the concept of cell type implicitly arises by the differential expression of factors and CAMs, which define the specific adhesive forces between cells. We derive the model by considering the forces acting on a population of cells with no proliferation or death to be conservative, which implies by mass conservation We assume that the cell dispersion velocity is proportional to the population density, which implies , (2.5) where is the dispersion constant.
The adhesion velocity vector depends on the adhesive bonds between the CAMs expressed in the cells and their neighbors within a sensing radius ( Figure 1a). This radius models the size of a cell, including their ability to reach and contact other cells through the cell body and through their protrusions such as filopodia. Following Newton's law and assuming that inertia is negligible for cell movements, the adhesion velocity vector is then inversely proportional to the cell radius (due to drag) and proportional to the sum of all adhesion forces between CAMs, such that where is the -dimensional unit spherical surface, is the radial distance, is the direction vector, , , , describes the nature of adhesive forces between CAMs and expressed from cell locations and , respectively, and their dependence on the cell density, and describes how the cell adhesive force depends on the radial distance . For simplicity, we assume 1 in this paper. The adhesive force between two CAMs expressed by two cells depends on their binding activity due to the CAMs relative concentrations within each cell. We assume that the binding activity follows the law of mass action, such as the adhesive force exerted on cells at location expressing CAM by cells at location expressing CAM depends on the product of their concentrations within their respective cells, given by where ℎ represents how the adhesive force depends on the local cell density. We assume a crowding capacity in the population limiting the cell movement due to adhesion towards dense regions, such as A nondimensional model is defined by rescaling with * , * , * , * , (2.10) and dropping the stars, we obtain the model

Simulations
We demonstrate the ability of the proposed model to recapitulate tissue shape behaviors due to the

Cell sorting behaviors
CAMs can bind to each other with different adhesive strengths, so cell-cell adhesion forces depend on the levels of expression of the different CAMs. These differences in cell-cell adhesion can cause an emerging cellular self-organization into different spatial patterns, a behavior demonstrated in vitro in a variety of animal cells, including amphibian (2), chick (35), zebrafish (18), and hydra (36).  The proposed continuous nondimensional model (2.11) can recapitulate these cellular sorting behaviors due to the differential expression of CAMs. Figure 3 shows four different sorting behaviors resulting from non-proliferating cells expressing either of two CAMs with different adhesive strengths (strength values as in (27)). All the simulations start with the same initial random configuration of disassociated tissue, where each initial aggregate contains cells expressing one of two different CAMs. Depending on the relative strength of self-and crossadhesion forces between the CAMs, the spatially-randomized tissues form aggregates that selforganize into patterns of engulfment, partial engulfment, mixing, or complete sorting. When the CAM self-adhesive strengths are asymmetric, similar to the retinal epithelial and retinal neural cells (Figure 2), the simulation recapitulates the engulfment behaviors observed in vitro ( Figure   3a). This sorting behavior is due to the differential expression of CAMs, where cellular aggregates expressing the CAM with stronger self-adhesive force (red) are tightly adhered, and hence surrounded by the cellular aggregates expressing the CAM with weaker self-adhesive force (green). However, when the self-adhesive strength of the two CAMs are equal, but still higher than the cross-adhesive strength, no cell aggregate is stronger than the other, and hence there is still sorting between the tissues expressing the different CAMs, but no engulfment (Figure 3b). In

Cell intercalation in proliferating cell populations
The spatial tissue behaviors in a population of proliferating cells can depend on both the expression levels and the adhesive properties of CAMs. This has been shown in in vitro assays of proliferating cell populations expressing similar or different nectin adhesion proteins (6). Figure 4 shows   We extend equations (2.2) to include a simple logistic cell growth term , such that where is the cell growth rate and is the cell carrying capacity. We assume that the daughter cells continue expressing the same CAMs than their parent cells, so that the relative CAM concentration in the daughter and parent cells are equal. In this way, we extend equation (

Dynamic regulation of adhesion during gastrulation
During gastrulation, Nodal acts as a diffusive morphogen forming a concentration gradient which induces mesendoderm differentiation (39). In zebrafish, Nodal is expressed in the yolk syncytial layer (YSL), a region of the yolk consisting of nuclei that have descended from the blastoderm (40). The YSL is divided into two segments: the internal YSL (iYSL), which is completely covered by the blastoderm, and the external YSL (eYSL), which protrudes beyond the blastoderm margin.
Only the nuclei of the eYSL are transcriptionally active, being the source of the Nodal signal that diffuses through a small area of the blastoderm at the region of the embryonic shield. All germ layers express similar levels of E-cadherin, but Nodal induces the up-regulation in expression of N-cadherin (41). These Nodal-induced cells with higher expression of N-cadherin have higher cell adhesion strengths compared to ectoderm cells (18) and they differentiate into mesendoderm.
Those not exposed to the Nodal gradient become ectoderm (42). Furthermore, the Nodal-induced cells with higher N-cadherin levels that differentiate into mesendoderm involute over the blastoderm margin towards the yolk (40). Figure 6 shows the zebrafish gastrulation between the 40% epiboly and shield stages, showing the involution movement of Nodal-induced cells towards the blastoderm margin. We extend the nondimensional model (2.11) to include a regulatory term for the CAMs and equations for the Nodal morphogen and other regulatory factors in the system, such that We employ this model to simulate zebrafish gastrulation due to the dynamic regulation of CAMs by the diffusion of morphogen Nodal (Figure 7). The model includes four CAMs: E-cadherin, Ncadherin, and those present in the EVL and yolk (labeled EVL and Yolk in Figure 7, for simplicity).
The adhesive constants between the CAMs are shown in Figure 7b. E-cadherin and N-cadherin have the same self-and cross-adhesion strengths, but E-cadherin has a stronger adhesion strength with the EVL in comparison to N-cadherin. The CAMs regulatory terms in , , are all zero, except for N-cadherin, which is regulated by the morphogen Nodal. Hence, the expression of N-cadherin depends on the levels of the morphogen Nodal, in addition to the cell density (more cells express more proteins) and a logistic saturation term, such that where is a regulatory constant, and is the concentration of Nodal. In addition to the morphogen Nodal, the model includes the regulatory factors expressed in the iYSL and eYSL regions of the yolk as two lumped variables in with zero diffusion constant and regulation.
Nodal diffuses and is expressed in the eYSL, hence its diffusion constant is not zero, and its reaction term depends on the eYSL factors in addition to natural degradation, such that , , , where is an expression constant, is the concentration of eYSL factors, and is the decay constant for Nodal.

Substituting (3.5) and (3.4) in (3.3), the zebrafish gastrulation model is defined with
(3.6)  including the expression and concentrations of CAMs. This dynamic regulation of CAMs is a key element in many biological processes, and to our knowledge has previously only been captured with hybrid models (25). However, numerically solving hybrid models are computationally infeasible for large cell numbers and their mathematical analysis is limited due to their discrete nature.
The proposed approach uses adhesion forces that arise dynamically from CAM expression, instead from explicit cell types. Single-cell force spectroscopy can directly measure the adhesion forces between cells expressing different levels of CAMs (45), which can be used to experimentally set the parameters of the model. For simplicity, the model assume that the binding forces are proportional to the product of the relative concentration of CAMs; however, more complex formulations of adhesion ligands and receptors are also possible (46). The extracellular matrix is an additional important element in cell-cell adhesion dynamics and it could be incorporated into the continuous model of adhesion (29). Together with adhesion forces, cell cortex elasticity and tension also can play a role in certain behaviors and be essential during tissue shape dynamics (18,47). Future work will extend the presented model to incorporate the role of these components in cellular behaviors.
The regulation of CAM expression and how these molecules affect large scale cellular behaviors is extraordinary important in both healthy and diseased states (48)(49)(50). The proposed model integrates genetic regulation of CAMs, the biophysical forces of adhesion that drive cell motion, and the subsequent cellular dynamics. This integrated view of dynamic adhesion will be essential for understanding tissue behaviors in developmental and cancer biology, as well as in bioengineering (51,52). The proposed continuum approach allows for the simulation of large cell populations (up to whole-organism scale), as well as the continuous phenomena involved in genetic regulation (e.g. morphogen gradients and CAM expression). Crucially, machine learning approaches for the reverse-engineering of the regulation of patterning (53-55) and cancer formation (56, 57) directly from formalized experimental data (58)(59)(60)(61) can be integrated with the proposed model, with the goal to discover the specific regulatory mechanisms of CAMs that give rise to key spatial phenotypes. In summary, the presented continuous modeling approach will pave the way for the understanding of the regulatory dynamics of cell adhesion essential in developmental, regenerative, and cancer biology.