The multi-levelled organizationof cell migration:from individual cells to tissuesbyHildur KnutsdottirB.Sc. Physics, University of Iceland, 2008M.Sc. Physics, Simon Fraser University, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Mathematics)The University of British Columbia(Vancouver)April 2018c© Hildur Knutsdottir, 2018AbstractCell migration is a complex interplay of biochemical and biophysical mechanisms. I investigatethe link between individual and collective cell behaviour using mathematical and computationalmodelling. Specifically, I study: (1) cell-cell interactions in a discrete framework with a spatialsensing range, (2) migration of a cluster of cells during zebrafish (Danio rerio) development, and(3) collective migration of cancer cells and their interactions with the extracellular-matrix (ECM).My 1D model (1), is approximated by a continuum equation and investigated using asymptoticapproximations, steady-state analysis, and linear stability analysis. Analysis and computations char-acterize regimes corresponding to cell clustering, and provide a link between micro and macro-scaleparameters. Results suggest that drift (i.e. due to chemotaxis), can disrupt the formation of cellularaggregates.In (2), I investigate spontaneous polarization of a cell-cluster (the posterior lateral line pri-mordium, PLLP) in zebrafish development. I use a cell-based computational framework (hybriddiscrete cell model, HyDiCell3D) coupled with differential equation model to track the segregationand migration of the PLLP. My model includes mutual inhibition between the diffusible growth fac-tors Wnt and FGF. I find that a non-uniform degradation of an extracellular chemokine (CXCL12a)and chemotaxis is essential for long-range cohesive migration. Results compare favourably withdata from the Chitnis lab (NIH).I continue using HyDiCell3D in (3) to elucidate mechanisms that facilitate cancer invasion. I fo-cus on: wound healing in a cell-sheet (2D epithelium), and cell-clusters (3D spheroids) embedded inECM with internal signalling mediated by podocalyxin, a trans-membrane molecule. Experimentaldata from the Roskelley lab (UBC) motivates the model derivation. I use the models to investigatethe role of cell-cell and cell-ECM adhesion in collective migration as well as the emergence of a dis-tinct phenotype (leader-cells) that guide the migration. ECM induced disruption in the localizationof podocalyxin on the cell membrane is captured in the model along with morphological changes ofspheroids. The model predicts that cell polarity and cell division axis influence the invasive poten-tial. Lastly, I developed quantitative methods for image analysis and automated tracking of cells ina densely packed environment to compare modelling results and biological data.iiLay SummaryI use mathematical and computational models to investigate coordinated cell migration (critical innormal development and cancer spread), focusing on the link between organization levels, fromindividual-cell to tissue-level behaviour, and from internal to external molecular networks. Startingfrom minimal models of interacting cells, I use analytical tools to elucidate the possible range ofmigratory behaviours. The analysis provides a link between individual and population level parame-ters. Furthermore, cell-cell mechanical and biochemical interactions coupled in a cell-based compu-tational framework facilitate investigation of the following systems: (1) For zebrafish sensory-organdevelopment, I model the chemical signalling involved in formation of distinct front and back cellswithin a cluster. (2) In cancer-cell sheets (2D) and clusters (3D) my models demonstrate that the mi-croenvironment can influence adhesion properties, polarity and cell division axis, leading to changesin cancer invasion. Lastly, I developed image analysis methods that enable quantitative comparisonof experiments and models.iiiPrefaceTheodore Kolokolnikov was a co-supervisor of the research presented in chapters 4 and 5.A version of chapter 6 has been published [Knutsdottir H, Zmurchok C, Bhaskar D, PalssonE, Dalle Nogare D, Chitnis AB, and Edelstein-Keshet L. (2017) Polarization and migration in thezebrafish posterior lateral line system. PLoS Comput Biol doi: 10.1371/journal.pcbi.1005451]. Iwas responsible for leading the team project, all major model formulations, 3D simulations andmanuscript preparation. Ajay Chitnis and Damian Dalle Nogare provided the biological data formodel design and validation. Dhananjay Bashkar and Eirikur Palsson were involved in the earlystages of the model formulation. Cole Zmurchok was involved in the model formulation, simu-lations of the 1D model and manuscript preparation. Leah Edelstein-Keshet was the supervisoryauthor involved in model formulation, model analysis and manuscript preparation.Dhananjay Bashkar improved the pre-processing of the experimental microscopy images for thesegmentation and cell tracking algorithm that I developed and used in chapter 7. Dhananjay alsomade scripts to automate the implementation of the tracking algorithm and data analysis for use inthe Roskelley lab.The experimental data in chapters 7 and 8 is from the Roskelley lab and the experiments wereconduct by Erin Bell, Pam Dean, Marcia Graves, Tak Poon and Alannah Wilson.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Research questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Research insights and relevance . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1 Computational background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1.1 Cell-based modelling frameworks . . . . . . . . . . . . . . . . . . . . . . 52.1.2 Computational framework comparison . . . . . . . . . . . . . . . . . . . . 62.2 Biological background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2.1 Collective migration in zebrafish development . . . . . . . . . . . . . . . . 72.2.2 Collective migration of cancer cells . . . . . . . . . . . . . . . . . . . . . 82.2.3 Coordinated migration of breast cancer cells and macrophages . . . . . . . 10v3 Computational model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123.1 Hybrid discrete cell model (HyDiCell3D) . . . . . . . . . . . . . . . . . . . . . . 123.2 Internal biochemistry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143.3 Chemical concentration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143.4 Cell orientation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143.5 Forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153.6 Cell motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153.7 Cell shape and proliferation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163.8 Model suitability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 Discrete model of migrating cells: analysis . . . . . . . . . . . . . . . . . . . . . . . 184.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184.2 Research questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214.3 Model derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214.4 Steady-state analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 234.4.1 Linear stability analysis of the discrete model . . . . . . . . . . . . . . . . 234.5 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 254.6 Continuum limit of the discrete model . . . . . . . . . . . . . . . . . . . . . . . . 254.6.1 Stability condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 274.6.2 Aggregation shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284.6.3 Numerical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 294.6.4 Slow speed regime: asymptotic approximation . . . . . . . . . . . . . . . 314.7 Alternative models and approaches for future work . . . . . . . . . . . . . . . . . 354.7.1 Symmetric model with directional bias . . . . . . . . . . . . . . . . . . . 354.7.2 Deriving the microscopic parameters from macroscopic observations . . . 364.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 375 Discrete model of interactions between distinct cell types: analysis . . . . . . . . . . 395.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 395.2 Research questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 395.3 Model derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 405.4 Steady-state analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 415.5 Continuum limit of the discrete model . . . . . . . . . . . . . . . . . . . . . . . . 445.5.1 Linear stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 455.6 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 465.6.1 The sensing radius and attraction between-species influence the dynamics . 465.6.2 Including a directional bias for one cell type . . . . . . . . . . . . . . . . . 485.6.3 Biological system: cancer cell-macrophage interactions . . . . . . . . . . . 48vi5.6.4 Hopf bifurcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 485.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 536 Collective migration in zebrafish development . . . . . . . . . . . . . . . . . . . . . 546.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 546.2 Research questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 556.3 Modelling frameworks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 566.3.1 Chemical signalling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 566.3.2 Cell-surface receptor levels . . . . . . . . . . . . . . . . . . . . . . . . . . 566.3.3 FGF and Wnt ligand concentration . . . . . . . . . . . . . . . . . . . . . . 596.3.4 Discrete cell model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 606.4 Tissue subdivision into leading and trailing zones . . . . . . . . . . . . . . . . . . 636.4.1 Theoretical considerations . . . . . . . . . . . . . . . . . . . . . . . . . . 646.4.2 The 1D spatial simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 666.4.3 Parameter dependence of the zone boundary . . . . . . . . . . . . . . . . . 676.5 Simulation results using HyDiCell3D . . . . . . . . . . . . . . . . . . . . . . . . 696.5.1 Formation of zones in 3D simulations . . . . . . . . . . . . . . . . . . . . 696.5.2 Migration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 716.5.3 Cell growth and proliferation . . . . . . . . . . . . . . . . . . . . . . . . . 746.6 Comparison with published experiments . . . . . . . . . . . . . . . . . . . . . . . 746.6.1 Laser ablation in silico . . . . . . . . . . . . . . . . . . . . . . . . . . . . 756.6.2 Comparison with Mutants and other experiments . . . . . . . . . . . . . . 766.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 787 The role of adhesion in collective migration of cell-sheets . . . . . . . . . . . . . . . 837.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 837.2 Research questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 847.3 Wound healing simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 857.4 Modelling results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 857.4.1 Changes in cell-cell adhesion, cell-ECM adhesion and the axis of cell divi-sion alter the leading-edge morphology of a cell-sheet . . . . . . . . . . . 867.4.2 Initial ratios of cell types affects the rate of motion of the leading edge . . . 877.4.3 Cell-cell adhesions affect relative distributions of cell types and leadingedge morphology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 887.5 Wound healing assay experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 897.5.1 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 917.6 Comparing wound healing simulations to experimental data . . . . . . . . . . . . . 917.6.1 Cell tracking algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93vii7.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 958 Coupling biophysical and biochemical processes to study cancer-cell spheroids . . . 988.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 988.2 Research questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 998.3 Biochemical model representing the podocalyxin dynamics . . . . . . . . . . . . . 998.3.1 A molecular switch model for podocalyxin . . . . . . . . . . . . . . . . . 998.3.2 RhoA dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1018.3.3 ODE model to describe podocalyxin localization within a cell . . . . . . . 1028.3.4 Sharp switch analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1038.4 Multi-cellular simulations coupling biochemistry to biophysics . . . . . . . . . . . 1068.4.1 HyDiCell3D simulation outline . . . . . . . . . . . . . . . . . . . . . . . 1068.4.2 Cluster morphology quantification . . . . . . . . . . . . . . . . . . . . . . 1078.5 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1088.5.1 The biochemical network influences cell polarity . . . . . . . . . . . . . . 1088.5.2 Variations in the total amount of podocalyxin in each cell alter cell polarity 1108.5.3 The micro-environment influences the cluster morphology . . . . . . . . . 1108.5.4 Axis of cell division plays a key role in cluster morphology . . . . . . . . . 1128.6 Comparing multi-cellular simulations with podocalyxin model to experimental results1138.6.1 Quantitatively comparing cluster morphologies . . . . . . . . . . . . . . . 1148.6.2 Comparing different extracellular matrices . . . . . . . . . . . . . . . . . 1178.6.3 Force generation of cells in contact with the extracellular matrix . . . . . . 1188.6.4 Manipulating the RhoA pathway . . . . . . . . . . . . . . . . . . . . . . . 1198.7 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1208.7.1 Modelling alignment and density of the ECM . . . . . . . . . . . . . . . . 1208.7.2 Cluster migration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1208.7.3 Impeding the RhoA activation experimentally . . . . . . . . . . . . . . . . 1218.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1229 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1259.1 Research summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1259.2 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1269.3 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128A Supporting Information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137A.1 The computational framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137viiiA.1.1 Chemical concentration . . . . . . . . . . . . . . . . . . . . . . . . . . . 137A.1.2 Forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137A.1.3 Cell shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138A.1.4 Cell division . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139A.2 Additional calculations for analysis in chapter 5 . . . . . . . . . . . . . . . . . . . 139A.3 Further details for zebrafish research . . . . . . . . . . . . . . . . . . . . . . . . . 141A.3.1 Model equations, derivation and scaling . . . . . . . . . . . . . . . . . . . 141A.3.2 Dose response experiment in zebrafish embryo . . . . . . . . . . . . . . . 143A.3.3 Parameter Estimation and Values . . . . . . . . . . . . . . . . . . . . . . . 143ixList of TablesTable 6.1 Parameters for model of collective migration during zebrafish development . . . 61Table 6.2 Parameters for 3D simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . 63Table 8.1 Cell cluster simulation parameters . . . . . . . . . . . . . . . . . . . . . . . . . 109xList of FiguresFigure 2.1 Signalling hypotheses in coordinated cancer cell migration . . . . . . . . . . . 11Figure 3.1 Computational framework flow chart . . . . . . . . . . . . . . . . . . . . . . 13Figure 3.2 Plot of interface force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16Figure 4.1 Signalling hypotheses in cancer migration . . . . . . . . . . . . . . . . . . . . 19Figure 4.2 Illustration of migration rules in discrete model . . . . . . . . . . . . . . . . . 22Figure 4.3 Instability condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24Figure 4.4 Single cell model simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 26Figure 4.5 Single cell model discrete and continuous comparison . . . . . . . . . . . . . 30Figure 4.6 Continuum equation for cells with directionality bias . . . . . . . . . . . . . . 30Figure 4.7 Comparison of analysis to simulations . . . . . . . . . . . . . . . . . . . . . . 34Figure 4.8 Continuum model simulations with directional signal . . . . . . . . . . . . . . 34Figure 4.9 Comparing analytic solution to simulations . . . . . . . . . . . . . . . . . . . 35Figure 5.1 Tumor cell, macrophage and endothelial cells signalling interactions . . . . . . 40Figure 5.2 Two species model (in)stability regions . . . . . . . . . . . . . . . . . . . . . 47Figure 5.3 Two species model (in)stability regions with drift . . . . . . . . . . . . . . . . 49Figure 5.4 Two species model comparing paracrine and within species attraction . . . . . 50Figure 5.5 Simulations with oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . 51Figure 5.6 Simulations of the chase-and-run mechanism . . . . . . . . . . . . . . . . . . 52Figure 6.1 The geometry and signaling of the PLLP . . . . . . . . . . . . . . . . . . . . 57Figure 6.2 Ligand-receptor dynamics forms a basis for the model. . . . . . . . . . . . . . 58Figure 6.3 Illustration of the cell-based simulations. . . . . . . . . . . . . . . . . . . . . 62Figure 6.4 Analysis of mutual inhibition dynamics leads to model insights. . . . . . . . . 65Figure 6.5 The Wnt-FGF mutual inhibition model . . . . . . . . . . . . . . . . . . . . . 67Figure 6.6 The position of the boundary between leading and trailing zones . . . . . . . . 68Figure 6.7 Experimentally observed Wnt and FGF . . . . . . . . . . . . . . . . . . . . . 70xiFigure 6.8 Simulations of boundary position between leading and trailing zones . . . . . . 71Figure 6.9 Successful migration in the discrete cell model . . . . . . . . . . . . . . . . . 73Figure 6.10 Abolished fibroblast growth factor (FGF) chemotaxis . . . . . . . . . . . . . . 74Figure 6.11 Effect of growth rate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75Figure 6.12 Recovery after partial ablation . . . . . . . . . . . . . . . . . . . . . . . . . . 76Figure 6.13 Additional laser ablation simulations . . . . . . . . . . . . . . . . . . . . . . . 77Figure 6.14 Dose response experiment using SU5402. . . . . . . . . . . . . . . . . . . . . 78Figure 7.1 Initial conditions in wound healing simulations . . . . . . . . . . . . . . . . . 86Figure 7.2 Invasion front for single cell type . . . . . . . . . . . . . . . . . . . . . . . . . 87Figure 7.3 Checkerboard pattern in the bulk of a cell-sheet . . . . . . . . . . . . . . . . . 88Figure 7.4 Various initial ratios between two distinct cell types . . . . . . . . . . . . . . . 89Figure 7.5 Wound closure time from simulations . . . . . . . . . . . . . . . . . . . . . . 89Figure 7.6 Result panel with various cell-cell and cell-ECM adhesion strengths . . . . . . 90Figure 7.7 Wound healing experiments with EpRas cells . . . . . . . . . . . . . . . . . . 92Figure 7.8 Tracks and displacement of simulated cells . . . . . . . . . . . . . . . . . . . 93Figure 7.9 Comparing experiments to simulations . . . . . . . . . . . . . . . . . . . . . . 95Figure 8.1 Single cell experiments using fibronectin bead . . . . . . . . . . . . . . . . . . 100Figure 8.2 Podocalyxin localization on the cell-membrane . . . . . . . . . . . . . . . . . 101Figure 8.3 Rho signalling pathway . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102Figure 8.4 Podocalyxin signalling pathway . . . . . . . . . . . . . . . . . . . . . . . . . 102Figure 8.5 Bifurcation diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104Figure 8.6 Sharp Switch Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105Figure 8.7 Time series from 3D simulations with Podocalyxin (podo) model . . . . . . . 107Figure 8.8 Illustration of biochemical and biophysical connection . . . . . . . . . . . . . 108Figure 8.9 3D simulations with different adhesion . . . . . . . . . . . . . . . . . . . . . 110Figure 8.10 Lumen formation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111Figure 8.11 3D simulations with various total amounts of podo . . . . . . . . . . . . . . . 112Figure 8.12 Morphology zones for different total amounts of podo . . . . . . . . . . . . . 113Figure 8.13 Simulations with and without the biochemical network . . . . . . . . . . . . . 114Figure 8.14 3D simulations with the podo model . . . . . . . . . . . . . . . . . . . . . . . 115Figure 8.15 Time series of podo simulations with proliferation . . . . . . . . . . . . . . . 116Figure 8.16 Phase image of experimental results . . . . . . . . . . . . . . . . . . . . . . . 117Figure 8.17 Segmented images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118Figure 8.18 Elongation index comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 119Figure 8.19 Additional simulations for ECM alignment . . . . . . . . . . . . . . . . . . . 120Figure 8.20 Additional simulations with cell tracks . . . . . . . . . . . . . . . . . . . . . . 121xiiFigure 8.21 Experimental results compared to simulations . . . . . . . . . . . . . . . . . . 123Figure A.1 Illustration of distances between cells . . . . . . . . . . . . . . . . . . . . . . 138xiiiGlossaryANOVA Analysis of Variance, a set of statistical techniques to identify sources of variabilitybetween groupsCAM cell adhesion moleculeCSF-1 colony stimulating factor 1, signalling moleculeCXCL12a stromal derived factor 1 (also known as SDF1a)CXCR4b chemokine receptor 4b that binds CXCL12aCXCR7b chemokine receptor 7b that binds CXCL12aE-cad epithelial cadherin a type of transmembrane protein that plays a role in cell adhesionECM extracellular matrixEGF epidermal growth factor, signalling moleculeEMT epithelial to mesenchymal transitionEpEN EpRas cells that express E-cad and N-cadEpE0 EpRas cells that express E-cadEp0N EpRas cells that express N-cadEpRas Pre-malignant mammary cells with epithelial characteristicsFGF fibroblast growth factorFGFR fibroblast growth factor receptorFN fibronection extracellular matrix protein that binds to integrin receptorsGFP green fluorescent protein used to label cells in experimentsxivHGF hepatocyte growth factorHyDiCell3D hybrid discrete cell-based modelling frameworkIC50 inhibition level required to reduce a function to 50% of its basal levelintegrin trans-membrane receptors that facilitate cell-extracellular matrix adhesionlef1 ranscription factor used as a marker for Wnt signalling domainmCherry red fluorescent proteins used to visualize cells in experimentsMDCK Madin-Darby canine kidney epithelial cell lineN-cad neural cadherin a type of transmembrane protein that plays a role in cell adhesionODE ordinary differential equationpea3 transcription factor used as a marker for FGF signalling domainPDE partial differential equationPLLP posterior lateral line primordiumpodo Podocalyxin, transmembrane anti-adhesive moleculeRhoA intracellular protein that plays an important role in cell polarity and contractilityWnt family of ligands that bind to the Frizzled receptorWntR receptors that bind to Wnt ligandxvAcknowledgmentsI sincerely thank my research supervisor Dr. Leah Edelstein-Keshet for her mentorship and endlessinspiration. I would also like to thank my committee members, Dr. Eric Cytrynbaum, Dr. JamesFeng and Dr. Calvin Roskelley, as well as the Keshet research group members for useful discussionsand debates. Special thanks go to my collaborators, Dr. Calvin Roskelley, Dr. Ajay Chitnis, Dr.Damian Dalle Nogare, Dr. Theodore Kolokolnikov and Dr. John Condeelis.Funding for my research was provided by the Canadian Cancer Society, the Natural Science andEngineering Research Council of Canada and the UBC Faculty of Science.Lastly, the support and encouragement of my family throughout this journey has been invalu-able. Thank you.xviChapter 1IntroductionDuring my lifetime I have a ∼50% chance of being diagnosed with cancer and a 25% chance ofthe disease killing me (based on the available statistics for Canada). The majority of cancer deathsare attributed to metastasis, the spreading of cancer through the body to distant organs. Frequently,metastasis is considered to mark a progression away from a curable disease because there are veryfew effective treatments available. Why is that? A tumor is a heterogeneous collection of cancercells as well as fibers, blood vessels and support cells [85]. Studying and understanding metastaticcancer requires methods that can capture interactions between these different factors. This hasproven challenging, leading to a lack of knowledge of fundamental mechanisms that underpin can-cer spread. One of these mechanisms involves interactions of cells that facilitate migration, whichis the main focus of my research.Cell migration is a complex interplay of biochemistry and biophysics, spanning multiple lengthand time scales. Cells can migrate individually or collectively, depending on their micro-environmentand their internal machinery (including genotype, molecular networks and adhesion molecule regu-lation). Environmentally influenced change in behaviour is referred to as phenotypic change: cellsthat start the same but undergo changes (such as morphological or biochemical), that result in dis-tinct cell behaviours. I will not attempt to treat the full complexity of cell migration but ratherfocus on interactions between cells and the coupling between intra- to inter-cellular processes. Forexample, in chapter 8, I couple a biochemical network inside each cell to adhesion to the so-calledextracellular matrix (ECM) in which the cells are embedded, to infer phenotypic changes on amulti-cellular level, observed experimentally.Mathematical and computational models can help elucidate key processes in cell migration byquantitatively linking individual cell behaviour to a multi-cellular system. Furthermore, modelsallow for a systematic comparison of different scales: from individual cells to tissues and fromexternal to internal pathways. The models are only as useful as the biological assumptions theyare based on are reasonable, and it is important to understand their limitations. Nevertheless, they1can be a valuable tool to distinguish between competing hypotheses, highlight knowledge-gaps ormissing mechanisms in a hypothesis, design experiments and make predictions.1.1 MethodologyMy research is interdisciplinary and conducted in collaboration with experimentalists. This ensuresa biological basis for model assumptions and establishes the interaction loop between modellers andexperimentalists that makes the results relevant. I start with a phenomenological description of cellmigration rules which allows for detailed mathematical analysis. I then extend my research to studyspecific biological systems. Although my initial motivation was to study collective cell migration inthe context of cancer metastasis, I found that much is still to be learned from well known migrationpatterns observed in development. Thus, my research spans processes from development to cancerand is organized from systems with the greatest availability of experimental data and hypothesesto onse where less is known. Moreover, the systems that I chose to study were contingent onavailability of laboratory data and suitable collaborators.Cellular level interactions are the main focus of my research, thus individual cells are the basicunit in my model derivation and design. The models are designed for a specific system to addressa set of open unanswered questions. Analytical tools can provide a fundamental understandingof regimes of behaviour for the cell migration properties. Next, I implement the interactions in amulti-cellular model to study patterns and morphologies of cell collectives on a macroscopic scale.An example of this method is in chapter 5 where I use linear stability analysis on a discrete modelfor cell-cell interactions between and within populations, and identify parameter regimes where cellclusters can form. Then, I numerically solve the model to show the different types of patterns thatcan emerge as well as the full range of dynamics.Furthermore, analytical results can become difficult to obtain as more of the complexity is in-corporated into the models, requiring the use of numerical methods. For example, in chapter 6, Iderive a model for how a cluster of cells, in zebrafish development, can segregate into front andback subpopulations with distinct fates. I use analytical tools, including phase planes, to map theregions. Then to study the migration of the cluster, additional levels of interactions are included ina numerical model. Specifically, I consider diffusible chemicals, forces between cells, and a binarydecision made by each cell subject to the state of its cell-surface receptors.Finally, I compare my simulation results to experimental data. In some cases (an example ofwhich is in chapter 7) this requires the development of image analysis and automated cell trackingmethods to extract quantitative data from biological experiments.21.2 Research questionsI study cell migration in systems ranging from development to cancer. The common themes can becaptured in the following overarching research questions:B How can observations of cell behaviour at one level of organization be used to infer propertiesat a different level of organization?B How does the micro-environment influence cell migration?In this thesis, I study both chemical and contact guided interactions between different cells andbetween cells and the ECM. Furthermore, I couple mechanical properties (including adhesion andpolarity) and biochemical pathways necessary for cell migration. Specifically, my research is guidedby the following more specific questions:• What can be inferred about individual cell properties (including proliferation and adhesion)from tissue level observations?• Given basic rules of cell-cell interactions, what can be inferred about tissue level behaviour?• How does polarity arise spontaneously and how is it maintained in a cluster?• How do changes in the interactions of cells with the micro-environment alter the morphologyof a cell cluster?• Are the proposed interactions with the ECM sufficient to explain migration and if not, whatcould be missing?1.3 Research insights and relevanceIn addition to gaining experience in mathematical modelling, numerical simulations and imageanalysis, in my research, I have learned some surprising and interesting facts about cell biology.Here, I share a few facts that have guided and motivated my research.a) There are few existing treatments that specifically target the metastatic potential of cancer.This may partly be due to the complex interactions that govern motility of cancer cells and thenotoriously difficult task of studying them both inside the body and in the lab. Therefore, there is asubstantial knowledge-gap when it comes to understanding the fundamental mechanisms involved.b) As revolutionary as genome sequencing is, it alone will not be delivering cures for all cancers.This is not to say that mapping various genes will not bring us closer to cures. Genomics can beused to identify driver mutations in certain types of cancer. However, cells with exactly the samegenome can be different phenotypically, which can result in drastic changes in behaviour. A majorfactor determining the cell phenotype is the micro-environment.3c) Cancer cells are extremely plastic. They can take up the characteristic of other types of cells.This sometimes prevents the immune system from detecting cancer cells or leads to drug therapyresistance.d) Lastly, I learned that developmental biology and cancer biology share many of the samesignalling pathways. Cancer cells can be difficult to study because they behave abnormally andundergo aberrant migration. In order to understand such abnormalities, it will be very useful tounderstand the fundamental mechanisms involved in the migration of normal cells.1.4 Outline In chapter 2, I will briefly introduce cell-based computational frameworks. I will also describethe biological systems that I study in more detail and introduce the relevant literature. In chapter 3, I will introduce a 3D cell-based computational framework that I use in partof my research. The framework simulates individual cells as deformable ellipsoids and theirmigration is based on forces. Internal biochemical signalling pathways can be captured insideeach cell and the concentration of external signalling molecules can be estimated by solvingreaction-diffusion equations on a grid. In chapter 4, I will design a 1D discrete model to study simple cell migration for a singlecell type subject to a directional drift. I will also derive an equivalent continuum model forcomparison and analysis. Chapter 5, will be based on the same framework as the previous chapter but I will apply it toa model of interactions between two cell types. In chapter 6, I will study collective migration during a specific process in normal developmentof zebrafish. I will derive a model for the polarization of a cluster of cells that leads topersistent migration over distances an order of magnitude greater then the cell cluster. In chapter 7, I will continue to study collective migration but with applications to cancercell invasion, focusing on cell-cell and cell-ECM adhesion properties. I will also investigatethe necessary mechanical attributes of leader cells, that is, cells with distinct properties thatappear to coordinate the collective migration of the follower cells. In chapter 8, I will focus on the role of adhesion and polarization on the shapes of 3D clustersof breast cancer cells. I will deduce the properties of individual cells that facilitate clusterdynamics, using a phenomenological model to infer the location of a specific trans-membranemolecule implemented in a multi-cellular system. Furthermore, interactions with the ECMwill influence cell polarity and proliferation leading to changes in cell cluster invasion. In chapter 9, I will summarize the research findings and outline future research directions.4Chapter 2Background2.1 Computational backgroundMany biophysical and biochemical cues act in unison to drive cell migration. These factors actacross multiple length and time scales and involve heterogenous cell populations. Models areuniquely capable of integrating the various interactions involved to provide insights, design ex-periments and suggest testable hypotheses that cannot be generated using experimental data alone,particularly when populations of cells are considered. Here, I briefly introduce cell-based simula-tion frameworks. The specific framework that I use to model cell migration will be described in thefollowing chapter.2.1.1 Cell-based modelling frameworksIn my research, I focus on cell-cell interactions and how individual cell behaviour influences thebehaviour of multiple cells that are either in physical contact or communicate through chemicalsignals. Modelling frameworks that use individual cells as the basic unit, are well suited to studythose interactions. Following is a short overview of the basic cell-based multi-cellular modellingframeworks. All of these frameworks can be used for hybrid simulations, meaning that they can becoupled with continuous models, for example for chemical concentrations or biochemical networks.Cellular Automata models are on-lattice models, meaning that space has been discretized intofixed lattices which the cell can occupy [25, 33]. The cell motion is usually based on a set of rulesinstead of physical properties. These simulations are fast and can incorporate large numbers ofcells. The lattice constrains the movement, cell-division and number of neighbours. Nevertheless,these simple models can often give reasonable results and they can sometimes be used to derivecontinuum equations for cell densities.Cellular Potts models are on-lattice models which originated in physics (the Potts model forinteracting spins on a crystalline lattice) [53, 84]. A cell is described as a collection of lattice sites5(e.g., squares or cubes) with shape changes governed by a Hamiltonian (energy minimization). TheHamiltonian can represent mechanical properties such as adhesion, and a Monte Carlo method isused to update the cell configuration. Calibrating a Monte Carlo Step (MCS) against actual time isa challenge and the underlying lattice can introduce discrepancies in computing the surface area (in3D) or perimeter (in 2D) of a cell.In Vertex based models the cell membrane is tracked using polygons represented with edges andnodes [28]. The movement of the vertices is determined based on mechanical forces. This permitsreasonable tracking of convex cell shapes, and actual forces between cells, but has the disadvantagesthat simulations have slower run-times and issues can arise during mesh restructuring (updating thelocations of edges and nodes while making sure that edges do not cross).Center-based models (CBMs) include off-lattice models with cells represented by their centerposition [21, 71, 77]. Some CBMs allow for variation or changes in cell shape and size, as wellas cell-cell interactions, cell-ECM interactions and proliferation. Under the same category, particlemodels where cell volume is neglected could be included. Cell-cell interactions are mainly definedbetween the centers of neighbouring cells. However, CBMs generally fail to capture the details ofcell shape.Publicly available cell-based simulation platforms include: NetLogo (Cellular Automata model)[98], Chaste (Cancer, Heart And Soft Tissue Environment which implements vertex model, cell-center based model and Cellular Potts model) [60], CompuCell3D (Cellular Potts model) [89] andMecaGen (Cell-based model that incorporates gene regulation) [19].2.1.2 Computational framework comparisonIn my research, I developed a type of Center-based model that is described in detail in the nextchapter. The framework was first developed by Eirikur Palsson to study slime mold Dictyosteliumdiscoideum migration patterns [70]. I have extended this framework to study cell migration andspecifically added interactions between cells and the microenvironment as well as feedback betweenbiochemical and biomechanical properties of cells. Briefly, the movement of cells is independent ofa lattice and governed by forces, this is different from Cellular Automata type of models. Detaileddescription of the cell boundary is not tracked in my model, contrary to Vertex based models. TheVertex based models can be problematic for tracking cell-cluster break-up or cell aggregation, bothof which are important in the systems that I will study. In my framework, the movement of cellsis determined by calculating forces that act on each cell. The forces can be compared directly tothe forces that are measured experimentally. Contrastingly, cell movement in Cellular Potts modelsis restricted to a lattice and governed by an energy minimizing function which requires carefulcalibration before it can be compared to experiments (to determine the real time equivalent of theaforementioned MCS).62.2 Biological backgroundIn this thesis, I study migration of cells in the developing zebrafish embryo and in the followingin vitro setups: 1D tracks, 2D cell-sheets and 3D spheroids. Following is a brief summary of therelevant biological background.2.2.1 Collective migration in zebrafish developmentThe zebrafish is a well studied model system. It has the advantage that the embryo is transparentmaking it easy to visualize. Because it is widely studied, many mutants and drugs that target specificfunctions during development are available. I will study a specific process in early developmentwhere a group of cells migrates from the head to the tail of the fish. This group of cells is calledposterior lateral line primordium (PLLP). During the migration, smaller clusters are left behind ata seemingly regular interval. These clusters later develop into the sensory hairs that the fish uses tosense the flow of water.The PLLP consists of more than 100 cells (≈4-5 cells wide) when migration is initiated, whichmigrate about 3 mm during the time interval 20-42 hours post-fertilization [34]. Cells in the trailingregion have epithelial characteristics and re-organize into rosettes which later become the sensoryhair cells. Generally, epithelial cells organize into cohesive sheets and attach to each other throughvarious junctions (adhesions) [39]. Contrarily, mesenchymal cells are characterized by multiple thinprotrusions and tend to migrate individually or be loosely cohesive. The cells in the leading regionare flat and have a mesenchymal morphology [49].Distinction between front and back in the cell clusterThe changes in cell morphology and behaviour across the PLLP is believed to be regulated throughthe Wnt (family of ligands that bind to the Frizzled receptor) and FGF (fibroblast growth factor)signalling networks [2, 49, 63]. The leading zone of the PLLP is characterized by high-levels ofreceptors that bind to Wnt ligand (WntR) activity while the trailing zone is dominated by fibroblastgrowth factor receptor (FGFR) activity. To the best of my knowledge, the specific family of ligandsthat bind to the Frizzled receptor (Wnt) ligand involved has not been identified although recentexperimental work by my collaborators (Damian Dalle Nogare and Ajay Chitnis at the NationalInstitute of Health) implicates Wnt10a. More research has been done on the downstream Wntsignalling network [2]. WntR activity at the front of the PLLP promotes the expression of FGF3and FGF10 ligands. It is also linked to increased expression of the gene sef which has been shownto inhibit FGFR activity. Consequently, the fibroblast growth factor (FGF) pathway is inhibitedin the leading cells, and the secreted FGF ligands can spread towards the rear. Activation of FGFsignalling in the trailing domain subsequently drives the expression of Dkk1, a diffusible inhibitorof Wnt signalling [2]. Therefore, the Wnt and FGF signalling domains are established in the leading7and trailing zones, respectively, of the PLLP.Directional migrationMigration is dependent on chemokine signalling across the PLLP and is guided by a strip of stromalderived factor 1 (also known as SDF1a) (CXCL12a), a chemokine that is produced by superficialmuscle pioneer cells [17, 38, 95]. Without CXCL12a the PLLP fails to migrate. However, theexpression level of CXCL12a is uniform and thus there must be another mechanism that providesthe directional cue for the PLLP. As it turns out, there is differential expression of two chemokinereceptors, CXCR4b (chemokine receptor 4b that binds CXCL12a) and CXCR7b (chemokine recep-tor 7b that binds CXCL12a), across the length of the PLLP which guides migration. Similar to theWnt-FGF polarity, the expression of chemokine receptor 4b that binds CXCL12a (CXCR4b) domi-nates in the leading end while expression of chemokine receptor 7b that binds CXCL12a (CXCR7b)dominates in the trailing zone with some overlap [17, 95]. The mechanisms that determine the dif-ferential expression of these chemokine receptors are connected to the development of the polarizedWnt-FGF signalling, but they are incompletely understood [2, 12].It has been proposed that cells expressing CXCR7b degrade CXCL12a at a faster rate than cellsexpressing CXCR4b. This sets up a local gradient of CXCL12a that ensures unidirectional mi-gration of the PLLP [88]. Additional complementary research suggests that the absolute level ofCXCL12a serves as the directional cue [16]. Moreover, recent laser ablation and cutting experi-ments in [16] showed that the trailing cells can chemotact towards a FGF ligand secreted by theleading PLLP cells. Evidence for this stems from several experiments described in [16]: (1) whentrailing cells are physically separated from the leading cells by laser ablation, trailing cells movetoward and eventually join with leading cells despite unpolarized expression of the chemokine re-ceptor CXCR4b in the trailing fragment. (2) This movement is inhibited by treatment with the FGFinhibitor SU5402, suggesting that FGF signaling is essential for this process. (3) Furthermore, anisolated trailing cell cluster generated by laser ablation can be induced to migrate directionally byplacing a bead soaked in recombinant FGF3 either ahead of or behind the cells. These experimentalresults suggest that while the leading and trailing cells coordinate migration by inducing a localCXCL12a gradient, the leader cells respond to CXCL12a with protrusive and migratory behaviourand the trailing cells undergo chemotaxis toward the leading cells [16].2.2.2 Collective migration of cancer cellsIt is commonly proposed that cancer progression towards a metastatic disease is initiated by epithelial-to-mesenchymal transition (EMT) where individual cells start breaking away from a primary tumor[30]. However, experimental examination of the tumor-stroma interface has led to observations ofcollective tumor buds. These buds consist of collectively migrating tumor cells which maintainsome cell-cell junctions and they correlate with metastatic disease. Direct evidence for collective8invasion of cancer cells was recently obtained by Cheung et al. [11]. Indirect evidence for collec-tive invasion is apparent from histopathological analysis of human cancers where clusters of cellshave crossed the tumor boundary and invaded the stroma without loosing cell-cell contact [37].When sheets (2D) or clusters (3D) of cells move collectively they have distinct properties, the cellsinteract and remain connected, multicellular polarity is maintained and the cells interact with themicroenvironment for example by structurally modifying the tissue or migration path [43].In my research, I focus on two specific experimental setups to study collective cancer cell mi-gration: wound healing assay (2D) and spheroid assay (3D).2D scratch wound healing experiments: emergence of leader cellsWound healing experiments are classical laboratory methods to study collective migration. Briefly,it involves plating cells in 2D and letting them organize to a densely packed cell-sheet, an epithe-lium. Next, various techniques can be used to make a wound in the epithelium: 1) removing cellswith a tool such a micro-pipette (scratch wound healing assay), 2) by removing a previously placedbarrier, or 3) removing cells using lasers, electricity or drugs [80]. Finally, the cells are observed,under a microscope, as they move to close the wound.At the invasion front of a migrating sheet or cluster of cells, leader cells that have distinctproperties from the rest of the cells can emerge [43]. Furthermore, leader cells are thought to bemajor drivers of collective invasion and thus they are worth characterizing. The leader cells aligncollagen fibres which form a path of least mechanical resistance for the cells [99]. They have alsobeen isolated experimentally and shown to have higher activity of certain molecules known to beinvolved in cell migration (Rac, β1-integrin and PI3K) [102]. The origin of leader cells remainsunclear and I will explore that in my research. Specifically, I will study whether cells with distinctadhesion properties can migrate to the leading edge of an epithelium and facilitate the migration.This study is motivated by experimental evidence showing that cells can migrate through a denseenvironment such as an epithelium [26].In wound healing assays, the invasion front has sometimes been shown to have a finger-likemorphology [76]. This morphology can change depending on the cell types and micro-environment.Similar fingering phenomena has been observed and modelled in other physical systems, i.e. Cahn-Hilliard equations [105].3D spheroids: podocalyxin localization and collective migrationPodocalyxin (podo) is a cell surface molecule with anti-adhesive qualities, formerly known asgp135. Podo was originally identified on cells in the kidney where it is essential for normal re-nal development, [64]. Using Madin-Darby canine kidney epithelial cell line (MDCK), it has beenshown that a single cell starts with podo located around the entire membrane [6]. When the celldivides the podo in the two cells is observed at the cell-ECM boundary but not at the cell-cell in-9terface. As a multi-cellular cluster is formed cell polarity is reoriented and the podo relocates tothe cell-cell interface (apical-basal polarity) to initiate lumen formation. The authors showed thatchanging the environmental condition (specifically by blocking the receptor β1-integrin) this re-location of podo gets disrupted. Other molecular perturbants were used to decipher the signallingpathway involved which I will elaborate on in chapter 8.My collaborators in the laboratory of Calvin Roskelley at UBC have shown that podo is an in-dependent predictor of poor prognosis in breast cancer patients [86]. The Roskelley group uses a3D experimental assay (spheroids) to study the invasive potential of breast cancer cells expressingvarious levels of podo. Spheroids are 3D cell clusters that are first grown on a 2D surface and thenembedded in different ECM, namely Matrigel (gelatinous protein mixture of basement membraneproteins that breast tumour cells must move through to invade the surrounding connective tissueswhere collagen-I is prevalent) or collagen-I [37]. In the Matrigel experiment, both control cell clus-ters and podo over-expressing cell clusters remain spherical. However, in the collagen-I experimentthe morphologies change and podo expressing cell clusters elongate significantly. The elongationof clusters can be an indicator of cells with invasive potential. I will study the morphology of theseclusters and analyze them with respect to a model that couples the podo relocation dynamics withthe cell mechanics.2.2.3 Coordinated migration of breast cancer cells and macrophagesIn order for cancer cells to spread to distant organs they need to first migrate into a blood vessel, aprocess known as intravasation. Research from the laboratory of John Condeelis at Albert EinsteinCollege of Medicine has shown that tumor cells and cells of the immune system, macrophages,coordinate their migration towards the blood vessels [15, 101]. These cells are not in physicalcontact (no adhesion) and their interactions are purely through chemical signals. In short, theseinteractions include a short ranged signalling loop involving epidermal growth factor (EGF) andcolony stimulating factor 1 (CSF-1). The tumor cells secrete CSF-1, which binds to CSF-1 receptorson the macrophages. In return, the macrophages secrete EGF, which binds to EGF receptors on thetumor cells. In addition, some types of tumor cells also have a receptor for CSF-1, thus they can beattracted to a signal from themselves. It has been shown that tumor cells interacting through bothEGF/CSF-1 and CSF-1/CSF-R signalling loops have an increased invasivity (in a study conductedin mice) [75]. Furthermore, recent results from the Condeelis lab suggest that hepatocyte growthfactor (HGF), a chemical secreted by cells that line the blood vessels, endothelial cells, acts as thedirectional cue for the tumor cells, [50]. These experiments were conducted in the lab using micro-patterned stripes of an adhesive substrate (fibronectin) which restrict the migration of the cancercells to 1D. Taken together, these experimental results suggest the signalling hypothesis illustratedin Fig. 2.1.10CSF-1EGFMacro- phageTumor cellBlood vesselHGFEndothelial cellsCSF-1EGFFigure 2.1: Proposed interactions between tumor cells, macrophages (cells of the immune sys-tem) and endothelial cells (cells that line the blood vessels). The endothelial cells pro-vide the directional cue towards the blood vessels by secreting HGF which the tumorcells can detect and migrate towards. The tumor cells are attracted to other tumor cellsvia a CSF-1/CSF-1-receptor signalling loop. Additionally, the tumor cells interact withmacrophages via an EGF/CSF-1 signalling loop, which enhances the directional migra-tion of tumor cells. EGF, CSF-1 and HGF are all diffusible signalling ligands.11Chapter 3Computational modelCell migration is a complex interplay between biochemical and biophysical attributes. It is furthercomplicated through interactions with other cells and the micro-environment. Choosing a modellingplatform depends on the spatial and temporal scales required to address the research questions. Themodel also needs to be able to capture the biological system in sufficient detail. The cell-basedmodelling platform that I use, described in this chapter, is classified as an off-lattice model. Off-lattice models can continuously track the cell location. Other cell-based models include CellularAutomata, Cellular Potts and Vertex based models described in the previous chapter.3.1 Hybrid discrete cell model (HyDiCell3D)The hybrid discrete cell model, or HyDiCell3D, is a cell-based model adapted from software (writ-ten in C) by Eirikur Palsson [70] with 3D output images and movies generated by openGL. Thecells are represented as deformable ellipsoids with finite volume. The interactions between cells areboth mechanical (i.e. cell-cell adhesion) and chemical (i.e. chemotaxis). Internal biochemistry ofeach cell is captured by ODEs describing the evolution of the factor and it can be linked to mechan-ical properties of the simulated cells. Random and chemotactic active forces govern cell motion.Attractive adhesive forces are associated with cell-cell and cell-substrate contact and a short-rangeexclusion force (sometimes referred to as repulsive force) prevents cells from overlapping. Themotion of cells is calculated from the net force acting on each cell.Fig. 3.1 illustrates a flow chart of the modelling framework. In the following sections, I will gointo the details of the model. In brief, once the initial conditions of the simulations have been set(such as number, types and location of cells), the following procedures are performed:1. Update internal biochemistry in each cell (ODEs for internal signalling pathways).2. Calculate all the extracellular chemical concentrations.3. Calculate the chemical concentration and gradients that each cell detects.12Calculate forcesUpdate chemical concentrationsMove all cellsOutput Make image filesChemotaxis, adhesionExclusionInitial conditions (Number of cells, cell types, properties, environment)CHAPTER 3. COMPUTATIONAL FRAMEWORK 15moving forward is slow while backwards sorting is faster andperhaps more efficient. This is consistent with my simulationssince backwards sorting of a few cells is like forward sortingof a large group of cells. The adhesive properties of Dd cellscould be manipulated by putting the expression of adhesionmolecules under the control of prestalk or prespore specificpromoters. If the prestalk cells are more adhesive than presporecells and the adhesion is not cell specific, the pattern shouldbe similar to Fig. 6(b), where prestalk cells clump together andmove anteriorly up the center of the slug. If the adhesion alsohas a cell specific component, then the two cell types willseparate more and form two distinct regions where the boundarybetween the prestalk and prespore cells is minimized, similar toFig. 5(b) or 6(b). If the prespore cells are more adhesive, moreprestalk cells would be found at the surface (Fig. 8 and 8(b)). Fromthe model I would predict that in mutants that do not express theprespore-specific adhesion molecules, the boundary between theprestalk and prespore zone would not be as well defined, with ahigher proportion of prespore cells in the prestalk zone than in thewild type.The sorting patterns that arise from the combined inter-actions of cell adhesion and chemotaxis, as shown here, arenot unique to Dd. There are many instances, such as duringgastrulation (Alberts et al., 1994), neural crest cell migration(Takeichi et al., 2000) and somitogenesis in zebrafish (Stickneyet al., 2000) where cells with differences in adhesive andchemotactic strengths sort out. The simulations of the modelcan thus be generalized and used to identify qualitative adhesiveproperties of cell types in those systems. The cell charac-teristics incorporated are general for most cell types but theparameter values may differ. Forces are carefully balanced in themodel, so it can be used to study the difference betweenmovement when cells get traction from the boundary versusneighboring cells.Since I was focusing on the cell sorting, I did not include theeffects of positional cell differentiation even though it clearlyplays an important role. Thus, improvements of the model wouldbe to include active cell differentiation, and to couple that withthe cell movements and signaling.AcknowledgmentsI would like to thank Leah Edelstein Keshet and Carol Wenzelfor useful suggestion during manuscript preparation. Supported inpart by NSERC Grant number RGPIN249789.Appendix AA.1. Cell propertiesThe cell is represented by its center of mass coordinates and bythe orientation and length of the principal ellipsoidal axes,denoted by the vectors a;b and c. The viscoelastic cells areassumed to be deformable ellipsoids with axes of length a, b and c.Each axis contains a nonlinear spring in parallel with a spring anda viscous element in series, also called a Kelvin model (Fung, 1993)(Fig. 11(a)), in line with findings of Chien (1984) and Skalak et al.(1984). k1, k2 are the spring constants and m1 is the viscosity of theelement. The fluid inside the cell, the cytosol, is essentiallyincompressible (bulk modulus 2! 109 N=m2), so the cell mustconserve volume under deformation. When one of the axis iscompressed, the other axes must be stretched, for volumeconservation, generating a counteracting force, represented as amodifying force, Fmod. The stress–strain relationship for the Kelvinmodel (Fung, 1993) including the Fmod gives us Eqs. (A.1) and (A.2)for the deformation of the i-axis of a cell for a given force on thataxis ðFiÞ.dridt¼ k1ðFi % FmodÞm1ðk1 þ k2Þþ dFi=dtðk1 þ k2Þ % rik1k2m1ðk1 þ k2Þ, (A.1)rarbrc ¼ ðra þ DraÞðrb þ DrbÞðrc þ DrcÞ ¼ Vellipse=ð43pÞ. (A.2)Here ri is the length of either the a, b or c axis in units of 10mm.Fmod is calculated from the volume constraint each timestep bysolving (A.1) to find Dra, Drb and Drc under the constraint of (A.2)(solution of a third-order polynomial for Fmod).A.2. ForcesA.2.1. Active forceI assume that when cell i actively tries to move, eitherrandomly or in response to a chemotactic signal, an active forceis generated. The cell attaches a ‘‘pseudopod’’ and applies thisforce to the substrate, or to that neighbor cell j in front of it whichhas the smallest angle between the vector rij (connecting the cellcenters of i and j) and the direction of the force. The cell alwaysorients its anterior–posterior axis (a-axis of the ellipsoid) in thedirection of the active force.When a cell attaches a pseudopod to a neighbor cell, the activeforce pulls the neighbor cell toward the first cell, while theopposing reaction force pulls the first cell forward. Thus, theactive force enters into the equation of motion for the second cell,and the reaction force enters into the equation of motion for theARTICLE IN PRESSκ1κ2µ1κ 1κ 2µ 1κ1κ1µ 1abrij j}idi }dj}drjriFactivebaFig. 11. (a) The combination of springs and dashpots (‘‘Standard Linear Solid Model’’) or the ‘‘Kelvin’’ model) that were used to model the viscoelastic properties of each ofthe ellipsoid axes. Here m1 is the viscous constant and k1 and k2 are the spring constants for Hookean springs. (b) Definition of the distance, d, between two cells. (c) Twooverlapping ellipsoids are used to represent close contact between neighboring cells. The cells are at equilibrium so the separation distance, d, is negative. The solid lineshows how real cells deform and the dotted lines indicate the ellipsoid shape.E. Palsson / Journal of Theoretical Biology 254 (2008) 1–1310Figure 3.1: Cells are repr sente s deformable ellips ids where each xis is r pr sented by aHookean spring with spring constant 1 in parallel with a Maxwell element consisting of viscouselement µ and a Hookean spring 2. Figure from Palsson [26].another one needs to be compressed in order to preserve the volume of the ell psoid. This generatesa modifying force, Fmod, which can be calculated by solving equations 3.1 and 3.2 simultaneously.dridt=1(Fi Fmod)µ(1 + 2)+dFi/dt(1 + 2) ri 12µ(1 + 2)(3.1)rarbrc = (ra +ra)(rb +rb)(rc +rc) = V/(43⇡) (3.2)where ri are the lengths of the diff ent axes of the ellipsoid, i can be a, b or c, Fi is the total f rceacting on axis i and ri is the change in length of axis i. When the force is removed the ellipsoidslowly relaxes back to its original shape.Before the cell begins to move it polarizes and establishes a front and back. In th model, toaccount for this the a axis of the ellipsoidal cells is always oriented in the direction the cell is movingin. If the concentration of a signaling molecule, either EGF or CSF-1, around a cell is above a setthreshold and the signaling molecule gradient is steep enough, the c ll will orient its a axis in thedirection of the gradient and start moving along it. The gradient calculated is the absolute gradientacross the cell diameter. If the gradient is below a set threshold the cells will chose a randomdirection with a biased toward the previous direction. When the un t-vector ~a of the ellipsoid isrotated in the new direction that the cell is moving in, the new ~b and ~c unit-vectors also need to becalculated. The new unit-vectors are calculated using the Gram Schmidt process for an orthogonalVisc elastic llsCHAPTER 3. COMPUTATIONAL FRAMEWORK 16Figure 3.2: Demonstration of how a cell (green sphere) can move independently of the lattice cubesu ed to c lculate c mical c ncentrations. Here, the green cell has volume i 8 diff rent latticecubes, which needs to be accounted for when calculating chemical concentration and secretion.basis where the old ~b and ~c vect rs are us d as th previous orthonormal basis. The lengths of thnew ~a, ~b, and~c vectors (ra, rb and rc) are calculated by finding where the new vectors would cut thesurface of the old ellipsoid. The force acting on a cell, further explained in section 3.4, determineshow much the ellipsoidal cells deform. The deformability of the ellipsoids helps visualize the cellspolarity and in practice represent the protrusion of cells along a chemical gradient. Using semi-hardspherical cells in my simulations would have produc d qu ntitatively similar r sults. However, sinceI hope to extend this work to simulate intravasation and extravasation of tumor cells, using ellipsoidalcells becomes important and does not add much to the computational time. For these simulations,the cells will need to squeeze through the walls of the blood vessels and the deformability of cellswill be necessary.3.2 Concentration of signaling moleculesThe tumor cells involved in the paracrine signaling loop with macrophages, secrete CSF-1 in re-sponse to a local EGF concentration above a set threshold. Similarly, macrophages secrete EGFin response to a high enough local CSF-1 concentration. The concentration of the two signalingC mical c entrations solved on a cartesian grid[A]: signalling molecule concentration. The secretion and degradation rates depend on the cells (N).Velocity depends on the viscosity (µ) of the mediumEquations of motionx: position vector v: velocity vectorGreen: epithelial cells Red: mesenchymal cellsGreen: tumor cells Red: macrophages Grey: fibronectionCells in a zebrafish embryo Yellow: FGF concentration Blue: Wnt concentrationBRIEF ARTICLE 3(30) Factive =8<: Fadh(x+ x0)e(x+x0)2) v0ex2i, if x > 0,Fcompress(x) 95 · ~rijk~rijk if x 0.@u@t= µ@2u@x2 @@x1(1 ✏(u+ v))u@v@x @@x2(1 ✏(u+ v))v@u@x(31)@v@t= µ@2v@x2 @@x3(1 ✏(u+ v))v@u@x @@x4(1 ✏(u+ v))u@v@x(32)u(x, t) = uss + ↵et+ikx(33)v(x, t) = vss + et+ikx(34)uss = vss = 0or(35)uss =1✏ v0(36)132ussvss > (µ 2uss)(µ 4vss)(37)where = 1 ✏(uss + vss)(38)(39)dudt= r · (D(1 + v)ruDurv u(1 u)(uara+ ubrb) + uv(vara+ vbrb)),(40)dvdt= r · (D(1 + u)rv Dvru v(1 v)(vara vbrb) + uv(uara+ ubrb)),(41)dadt= r2a+ u a,(42)dbdt= r2b+ v b.(43)@[A]@t=Diffusionz }| {Dr2[A] +Secretionz}|{Psec Degradationz }| {(kglobal + klocal)[A](44)dLdt= kL(45)dCdt= vrC(46)dFdt= Dr2F F + ↵C2 THE AUTHORFadhesive = Fadh(x+ x0)e(x+x0)2) v0ex2i, x 0(16)Frepulsive = Fcomp(x)3/2, x < 0(17) =rcell2✓1di+1dj◆(18)CSF-1 concentration:d[A]dt=Diffusionz }| {Dr2[A] +Secretionz}|{lsec Degradationz }| {kdeg[A]EGF concentration:dEdt= Dr2E + Esec(C) kegfE(19)(20) Csec =Xp2(ijk)SijkcellpScellkt morsec H(Ecellp Eth)Ecell p1 + Ecellp(21) Esec =Xq2(ijk)SijkcellqScellkmacrosec H(Ccellq Cth)Ccellq1 + Ccellqkcsf = kext +0@ Xq2(ijk)SijkcellqScell1A kcsflocal,(22)kegf = kext +0@ Xp2(ijk)SijkcellpScell1A kegflocal.(23)(24) Factive =⇢Fchemotax if gradient steepness is above threshold,Frandom otherwise.(25) Factive =⇢Fchemotax above threshold,Frandom otherwise.(26) Fcell = Factive +XN(Fadhesion + Frepulsive)(27) Frepulsive =(0 if x > 0,Fcompress(x) 95 · ~rijk~rijk if x 0.dxidt= vi,(28)vi =Fcellµ.(29)2 THE AUTHORFadhesive Fadh(x+ x0)e(x+x0)2) v0ex2i, x 0(16)Frepulsive = Fcomp(x)3/2, x < 0(17) =rcell2✓1di+1dj◆(18)CSF-1 concent ation:d[A]dt=Diffusionz }| {Dr2[A] +Secretionz}|{lsec Degradationz }| {kdeg[A]EGF concentration:dEdt= Dr2E + Esec(C) kegfE(19)(20) Csec =Xp2(ijk)SijkcellpScellktumorsec H(Ecellp Eth)Ecellp1 + Ecellp(21) Esec =Xq2(ijk)SijkcellqScellkmacrosec H(Ccellq Cth)Ccellq1 + Ccellqkcsf = kext +0@ Xq2(ijk)SijkcellqScell1A kcsflocal,(22)kegf = kext +0@ Xp2(ijk)SijkcellpScell1A kegflocal.(23)(24) Factive =⇢Fchemotax if gradient steepness is above threshold,Frandom otherwise.(25) Factive =⇢Fchemotax above threshold,Frandom otherwise.(26) Fcell = Factive +XN(Fadhesion + Fexclusion)(27) Frepulsive =(0 if x > 0,Fcompress(x) 95 · ~rijk~rijk if x 0.dxidt= vi,(28)vi =Fcellµ.92 THE AUTHORFadhesive = Fadh(x+ x0)e(x+x0)2) v0ex2i, x 0(16)Frepulsive = Fcomp(x)3/2, x < 0(17) =rcell2✓1di+1dj◆(18)CSF-1 concentration:d[A]dt=Diffusionz }| {Dr2[A] +Secretionz}|{lsec Degradationz }| {kdeg[A]EGF concentration:dEdt= Dr2E + Esec(C) kegfE(19)(20) Csec =Xp2(ijk)SijkcellpScellktumorsec H(Ecellp Eth)Ecellp1 + Ecellp(21) Esec =Xq2(ijk)SijkcellqScellkmacrosec H(Ccellq Cth)Ccellq1 + Ccell qkcsf = kext +0@ Xq2(ijk)SijkcellqScell1A kcsflocal,(22)kegf = kext +0@ Xp2(ijk)SijkcellpScell1A kegflocal.(23)(24) Factive =⇢Fchemotax if gradient steepness s abov thre hold,Frandom otherwise.(25) Factive =⇢Fchemotax above threshold,Frandom otherwise.(26) Fcell = Factive +XN(Fadhesion + Fexclusion)(27) Frepulsive =(0 if x > 0,Fcompress(x) 95 · ~rijk~rijk if x 0.dxidt= vi,(28)vi =Fcellµ.(29)Update internal biochemistryBRIEF ARTICLE 3(30) Factive =8<: Fadh(x+ x0)e(x+x0)2) v0ex2i, if x > 0,Fcompress(x) 95 · ~rijk~rijk if x 0.@u@t= µ@2u@x2 @@x1(1 ✏(u+ v))u@v@x @@x2(1 ✏(u+ v))v@u@x(31)@v@t= µ@2v@x2 @@x3(1 ✏(u+ v))v@u@x @@x4(1 ✏(u+ v))u@v@x(32)u(x, t) = uss + ↵et+ikx(33)v(x, t) = vss + et+ikx(34)uss = vss = 0o(35)uss =1✏ v0(36)132ussvss > (µ 2uss)(µ 4vss)(37)where = 1 ✏(uss + vss)(38)(39)dudt= r · (D(1 + v)ruDurv u(1 u)(uara+ ubrb) + uv( vara+ vbrb)),(40)dvdt= r · (D(1 + u)rv Dvru v(1 v)(vara vbrb) + uv(uara+ ubrb)),(41)dadt= r2a+ u a,(42)dbdt= r2b+ v b.(43)@[A]@t=Diffusionz }| {Dr2[A] +Secretionz}|{Psec Degradationz }| {(kglobal + klocal)[A](44)dBdt= f(B, ~y) g(B, ~y)(45)dLdt= kL(46)dCdt= vrCB: describes a biochemical state of cell. y: can be other internal biochemical factors, or from the environmentFigure 3.1: Illustration of the computational framework. After the initial conditions have beenset to represent a specific biological system or experiment the biochemical and mechan-ical properties are calculated for each cell. At the end of each time step all the cells aremoved to their new location and the information is saved to a file. At the end of the sim-ulations, OpenGL is used to generate image and video files to illustrate the simulationresults.134. Orient cells either based on the chemical gradient that it detects or the internal biochemistry,depending on the system that is being modelled.5. Calculate the forces acting on each cell.6. Update all cell locations.7. Repeat these steps.3.2 Internal biochemistryThe model can track internal biochemical factors such as amount of a protein, gene expression orthe proportion of activated receptors. In each cell, an ODE is solved that describes the state of someinternal biochemical factors, B,dBdt= f (B,y)−g(B,y), (3.1)where y represents a vector of other internal or external factors that influence the dynamics of B.External factors that affect B can include mechanical forces, the biochemical state of a neighbouringcell or interactions with the environment. Similarly, the biochemical factor B can influence themechanical properties of the cell i.e. the strength of cell-cell adhesion can be determined based onthe value of B.3.3 Chemical concentrationTo track ligand concentrations, reaction-diffusion PDEs are solved on a 3D cartesian grid (witheach grid cube representing (10 µm)3) for each ligand. The cells are not confined to move on thisgrid (off-lattice migration of cells) however they can act as sources and sinks for external ligands(signalling chemicals). This is accounted for in the model by calculating the fraction of the cell’ssurface area located in each lattice cube and distributing or depleting the ligand amount accordingly.The cells can react to the local gradient concentration of an external ligand. Thus, the finitesensing size of cells (normally assumed to be the cell size), that can possibly overlap several voxelsin the chemical grid, must be accounted for. I take an average over those voxels, to calculate thedetected concentration and use the polarity axis to calculate a gradient by linearly interpolatingbetween the concentration in front and behind the cell.3.4 Cell orientationFor chemotacting cells the major axis of the ellipsoid, the a-axis, is oriented up the chemical gradientto capture the cell polarity. In the absence of a chemotactic signal, the cells orient in a randomdirection but a bias can be specified for moving towards the direction they were previously orientedin. Depending on the system being modelled, factors other than chemotactic gradients can influence14the polarity of the cells. For example, the polarity of neighbouring cells, an internal biochemicalfactor or elements in the environment, such as structural fibers.3.5 ForcesTo calculate the movement of the cells, I first need to calculate the forces acting on each cell. Thefollowing is a summary of the forces on and between cells.Active force, Factive, represents the traction force from the environment. It can be due to chemotaxisgradients, adhesion to environmental components, gradients of substrates rigidity or spontaneous(random) active propulsion. The chemotaxis force has a fixed constant magnitude whenever thelevels of a ligand are detectable (i.e. above a specified threshold). The direction of the active forceis determined from a vector sum of the local ligand gradient. Failing to detect a chemical gradient,the cell applies a weaker force in a random direction but biased towards the direction from theprevious time step.Substrate adhesion force, is only in the z direction, acting to attach cells to a 2D substrate; thisforce is orthogonal to the xy plane and it has no effect on migration of cells in the xy plane, nor ondrag forces (which act to oppose that migration).Interface forces, include both adhesion and volume-exclusion (attraction and repulsion) at cell-cellinterfaces. Since ellipsoids are not space-filling, they are allowed to overlap somewhat at equilib-rium to capture densely backed sheets or clusters of cells. This overlap also allows neighbouringcells to have a finite contact area, which would not be the case for impermeable rigid ellipsoids.Further overlap is avoided by a power-law exclusion force that depends on x, the distance betweenthe cell surfaces along the vector connecting their centers. The adhesion force between two cellsacts along the same axis, with a magnitude that goes to zero for cells that overlap or that are farapart. This force pulls cells towards their equilibrium spacing, but acts locally. The interface force,Finter f ace, is plotted in Fig. 3.2.3.6 Cell motionThe motion of each cell is determined by the net force Fcell , defined as:Fcell = Factive+∑NFinter f ace, (3.2)where N is the number of neighbours. The position of the cell is updated by the (low Reynold’snumber regime) equation of motion,d~xdt=Fcellµ, (3.3)where µ is a drag coefficient.15Figure 3.2: An illustration of the shape of the interface force, Finter f ace, as a function of thedistance between cells (see Eqn. (A.2)).3.7 Cell shape and proliferationThe cell is modelled as a viscoelastic ellipsoid with semi-minor axes a,b and c. The viscoelasticity isobtained by representing each axis with the Standard linear solid model which consists of a Maxwellelement (a spring, with spring constant κ1, in series with a dashpot whose frictional coefficient is µ1)in parallel with a spring (with spring constant κ2), illustrated in Supporting Information Fig. A.1.(See [70] for typical mechanical parameters.) The forces acting on a cell are resolved onto thethree ellipsoidal axes and the shape of the cell is obtained from these mechanical elements undervolume conservation. The deformation capabilities of the cells can capture the commonly observedelongation of migrating cells. Additionally, the elongation of cells alters the cell-cell interactionsbecause the adhesion force between cells depends on their contact area. However, real cells are notellipsoidal and this framework cannot capture complex protrusive cell shapes.The cells are able to undergo division (mitosis) in the simulations. The cells grow in volumeover time with a constant growth rate for each of the three primary axes of the ellipsoids. Whena cell has doubled its volume it divides. The axis of cell division can be chosen at random or theelongated cell-shape can be used to determine division axis. The two daughter cells generally inheritthe parameter values from the mother cell. See Supporting Information A for more details on cellgrowth.163.8 Model suitabilityThis modelling framework has the advantages of being fast, simulating thousands of cells in minuteson a laptop, fully 3D and having the flexibility to design simulations to capture different biologi-cal experiments. The cells are deformable to capture overall aspects of migrating cell shapes andvariability in the number of neighbours. Therefore, this framework is well suited when studyingdensely packed cells where the detailed shape of individual cells is not the focus of the research butrather the collective behaviour, such as migration, polarity and tissue morphology. Cell migration isa very complex process which incorporates events at different size and time scales. This model caninclude interactions between cells that happen on the mesoscopic scale and couple them to internalbiochemical events that happen on the microscopic scale. The HyDiCell3D is also a good choice formodelling migration of multiple cells since the cells are discrete (off-lattice model) and their move-ment is based on calculating forces acting on each cell. An added benefit of using a force-basedmodel is that it can be parametrized using biophysical and biochemical parameters measured exper-imentally. Results from cell-based models, like this one, can be compared directly to experimentalobservations.17Chapter 4Discrete model of migrating cells:analysis4.1 IntroductionThe main objective of my research is to understand interactions of cells with their micro-environmentin the context of cell migration. These include interactions with different cell types within closeproximity. Mathematical models can sometimes help advance the fundamental understanding ofcell-cell interactions. In this chapter, I derive a simple one-dimensional (1D) model that can bestudied both analytically and numerically. In some cases, 1D models are a good approximation ofthe biological system [4, 41, 61]. Experimental setups, such as micro-fluidic devices and micro-patterned substrates, restrict the migration of cells to 1D. Furthermore, two dimensional (2D) ex-periments in well-mixed conditions can be represented in 1D by averaging the cell density over oneof the spatial dimensions.In terms of cancer biology, experiments have shown that cancer cells and cells of the immunesystem, macrophages, migrate together along collagen fibers towards blood vessels [35, 50, 100].The restriction of the cells to move along these fibers makes exploring 1D mathematical models ofmigration a good approximation, see Fig. 4.1. All the known and hypothesized interactions of thesecell types have not been accounted for in this chapter. However, I will assemble a framework, bothin this chapter and the following one, to study these interactions and I explore the insights that theframework can provide. Specifically, I will design a framework where the parameters for individualcell behaviour can be linked to those of a system of cells migrating collectively.Mathematical models for cells come in many varieties including discrete and continuum mod-els as well as hybrid models which combine both approaches. In discrete models, cell motion isdescribed based on a set of rules for each cell. The rules are made to capture microscopic behaviourof individual cells. The model parameters and results can easily be compared to results from exper-18Tumor cellBlood vesselHGFFigure 4.1: Proposed interactions between tumor cells and endothelial cells (cells that linethe blood vessels). The endothelial cells provide the directional cue towards the bloodvessels by secreting HGF which the tumor cells can detect and migrate towards. Addi-tionally, the tumor cells are attracted to one another through a chemical signalling loop(CSF-1/CSF-1R).iments. However, analyzing parameter dependence in these models is not always straight forwardand may require numerical methods. In continuum models, individual behaviour of cells is notdescribed but a well mixed assumption is made and the macroscopic behaviour of many cells ismodelled. This modelling approach is good when studying dynamics of a large number of cellsand multiple analytic techniques have been developed for these models [41]. In principle, these twomodelling approaches are comparable. However, the connection is not straight forward and it posessome challenges, such as: How can a probability of lattice occupancy be compared to a local den-sity? What are the characteristic spatial and temporal scales? How many terms in a Taylor seriesshould be kept when expanding model equations? What limits in the discretization of space andtime (∆x and ∆t) should be taken? [61, 66]Here, I derive a discrete model for a specific type of cell movement and obtain an equivalentcontinuum model for comparison. By deriving a continuum limit of the discrete model, analyticaltools can be used to study the model dynamics and the parameter sensitivity. This research isbased on an existing framework introduced by Chavy-Waddy and Kolokolnikov, to study bacteriamovement [9]. I will explore a variant of this model that includes a directionality bias and apply itto cell-cell interactions.The discrete model in [9] was based on experimental observations as well as previous modellingwork that derived ODEs for the migration of left and right moving bacteria from local interactionsdescribed in a probabilistic framework [31].In the model in [9], U j represents the density of bacteria in a given bin, j, and Wj is the sumof the densities of right movers in bin j+1 and left movers in bin j−1. (Here I have changed the19notation from [9], Wj = Vj, to avoid confusing the variable V with the rate v of biased directionalmigration for the cells in my model.) In each time step, cells either keep moving in the samedirection (with rate a) or they switch to a direction of higher cell density (and move at a rate c). Inthe mean-field limit for the bacteria model, the resulting ODEs (Eqn. (8) in [9]) were:dU jdt= a(U j−1+U j+1−Wj)− (a+ c)U j + c(U j−1η+j−1+U j+1η−j+1), (4.1a)dWjdt= (a+ c)(U j−Wj), (4.1b)where η±j =∑dk=1U j±k∑dk=1(U j+k +U j−k). (4.1c)d is the sensing distance in number of bins (unit-less parameter) and η±j is a nonlinear term thatdescribes the fraction of all bacteria within d bins of j that are to the left (-) or right (+). Definingη±j in this way leads to η+j +η−j = 1. Using linear stability analysis of the homogeneous steadystate, Chavy-Waddy and Kolokolnikov found the following threshold for the onset of instability:c0 =2ad. (4.2)When the rate c> c0 aggregation occurs in the system. This result suggests that larger sensing radiusaids the onset of aggregation. Conversely, highly persistent motion (large a) makes aggregationmore difficult. Next, the authors assumed that the inter-bin distance is h and let U j(t) = u(x, t) sothat U j+k(t) = u(x+kh, t) and similarly Wj(t) =w(x, t). By taking the limit of large number of bins,nondimensionalizing the equations and setting u(x, t) = w(x, t) (because they found that, except fora brief transition period, w quickly approaches u) the authors derived the following continuum limitof this model:ut =−uxx−uxxxx+α(uxuxxu)x(4.3)where α is an aggregate parameter that depends on all the parameters of the discrete model, a,c andd in the following manner:α =c(2d+1)(d+1)2c(1+d(d2+2d+3))−2a . (4.4)By solving the steady state of the PDE in Eqn. (4.3), the profiles of the aggregates could be calcu-lated, namely:u(x, t) =C[sech(√α−12x)] 2α−1. (4.5)The authors found good agreement between numerical and analytical results. The work in [9]20provides a novel framework to compare discrete and continuum models.In order to derive a discrete model for cell-cell interactions, I will first need to address thefollowing: How should cell-cell interactions be represented? How can I infer reasonable rules fora model from observed cell behaviours? What features can and should be included in the model?Similarly, which features should be left out for now? What behaviour am I interested in capturing?In this chapter, I will present an extension of the framework from [9], where I have includeda directionality bias (or drift term) to capture a type of directional migration of the cells due toenvironmental conditions. I will derive a continuum limit and study how the model parametersaffect the dynamics of the cells.4.2 Research questionsHere, I extend the framework from Chavy-Waddy and Kolokolnikov, [9], to study cell-cell interac-tions. The research questions that I explore in this chapter include:1. Under what parameters do the cells form aggregates?2. How does directional bias (i.e. due to chemotaxis or haptotaxis) affect aggregation of cells?3. How do cell-cell interactions affect aggregation?4. Can I predict the rate at which aggregates move?5. Can I derive a continuum limit for the discrete model? If so, how do the parameters of thediscrete model compare to the parameters of the continuum model?4.3 Model derivationI start by deriving a discrete model of the density for one cell type, U , with a directional bias of ratev, to the left. This bias can be interpreted as a bias in cell migration due to environmental conditionssuch as chemotaxis (gradient of a chemical) or haptotaxis (gradient in adhesion binding sites). Idescribe the movement of cells that are restricted to move in 1D along a line that has been splitinto intervals or bins of a finite size ∆x. The index j = 1, ...,n refers to the bin x j = j∆x. In thisformalism a typical bin size would be on the order of a cell diameter, ∆x= 5−20µm. For individualcell migration of speed 1µm/min, a characteristic time step for ∆x = 10µm should be 10 minutesso that cells typically make transitions that are no larger than one bin size per time-step. Goingforward, I will use a non-dimensional form of the equations to keep this phenomenological modelof cell migration general, that is, I rescale time by the above characteristic time step, and distanceby the typical distance that a cell moves during that characteristic time.A cell in bin j can move to bins j− 1, j+ 1 or stay in bin j with probabilities that depend onthe local cell density. The probabilities need to be small to ensure that a cell does not move more21than one bin size in a given time step and they can depend on the size of the bins. In my model, thedirectional bias, v, acts to increase the probabilities of moving to the left. The following is assumedabout the dimensionless transition rates: a is the rate at which the cell moves randomly, c is therate at which the cell moves after choosing a direction based on neighbouring cell densities and theprobability of a density dependent cell directionality is defined in Eqn. (4.1c). An illustration oftransitions between bins can be seen in Fig. 4.2. These assumptions lead to the balance equation:dU jdt=U j−1(a+ cη+j−1)+U j(−2a− c− v)+U j+1(a+ v+ cη−j+1). (4.6)The difference between this model and the previous model in Eqn. (4.1) is (1) the directionalbias term v and (2) since this model was not derived from left moving and right moving cells Ihave assumed U =W from Eqn. (4.1) (as was done in [9] when the continuum limit was derived,Eqn. (4.3)). Note that because of the way I defined and implemented the directional bias to the left,the balance equation is asymmetric.! ! ! ! ! ! ! ! ! ! ! ! !j-1 j j+1︸Δxd=1vBRIEF ARTICLETHE AUTHOR(1)dnidt= koffni + koni(na ni)(2) T+j1 = a+c4fy(x) = Q(y, x)Q(y, y)(3)= R(x) + g(x) b(x) + s g(x)(4)= R(x) b(x) + s(5)R(x) =x1 + x2(6)b(x) = x(7)(8) D(x) =@fy(x)@x=@R(x)@x @b(x)@xEvolutionary stable: D(x⇤) = 0Convergent stable: dDdx |x=x⇤ < 0@I@t= ↵(G)P 1I + 1(G)I @@x(I@G@x)(9)@P@t= ↵(G)P 2P + 2(G)P(10)@G@t= D@2G@x2+ s1(I)I 3(I + P )(11)@@t= µ@2@x2 @@x(1@c@x)(12)@T@t= µ@2T@x2 @@x(2T@e@x)(13)@c@t= D@2c@x2+ s1T 1c(14)@e@t= D@2e@x2+ s2 2e(15)1BRIEF ARTICLETHE AUTHOR(1)dnidt= koffni + koni(na ni)(2) T+j = a+c4fy(x) = Q(y, x)Q(y, y)(3)= R(x) + g(x) b(x) + s g(x)(4)= R(x) b(x) + s(5)R(x) =x1 + x2(6)b(x) = x(7)(8) D(x) =@fy(x)@x=@R(x)@x @b(x)@xEvolutionary stable: D(x⇤) = 0Convergent stable: dDdx |x=x⇤ < 0@I@t= ↵(G)P 1I + 1(G)I @@x(I@G@x)(9)@P@t= ↵(G)P 2P + 2(G)P(10)@G@t= D@2G@x2+ s1(I)I 3(I + P )(11)@@t= µ@2@x2 @@x(1@c@x)(12)@T@t= µ@2T@x2 @@x(2T@e@x)(13)@c@t= D@2c@x2+ s1T 1c(14)@e@t= D@2e@x2+ s2 2e(15)1Tj+1 = a+c2+ ⌫Tj = a+3c4⌫Figure 4.2: Illustration of the transitions between bins for a special case to demonstrate themigration rules in the discrete model. Each cell is represented with a red circle. Thetransitions between bins are represented with T . In this example, d = 1, η+j =14 andη−j =34 .I set out to address the following questions:• In this simplified model, will clusters emerge? If so, what are the necessary conditions for theonset of clustering?• Can the shape of the clusters be predicted?• What determines the rate at which aggregates move?224.4 Steady-state analysisThe steady states of Eqn. (4.6) are found by setting the left hand side equal to zero. The system hasa homogeneous steady state, U =U0, where U0 is a constant. Negative values for cell density do nothave a biological application, thus U0 ≥ 0. For homogeneous steady states, the η±–terms are equalto one half.To understand the dynamics of this model, I start by analyzing the stability of the steady states.The interesting dynamics will occur outside the stable steady states.4.4.1 Linear stability analysis of the discrete modelTo study the stability of the homogeneous steady states, U0, of Eqn. (4.6), I introduce a smallperturbation around U0, namelyU j =U0+φ j,where |φ j| U0. I substitute this perturbation into Eqn. (4.6), which leads todφ jdt= (U0+φ j−1)(a+ cη+j−1)+(U0+φ j)(−2a− c− v)+(U0+φ j+1)(a+ v+ cη−j+1).(4.7)I then use a Taylor expansion and only keep linear terms in φ j/U0 to make the following approxi-mation:η±j =∑dk=1U0+φ j±k∑dk=1(2U0+φ j+k +φ j−k)≈(12+∑dk=1 φ j±k2dU0)(1− ∑dk=1 φ j+k +φ j−k2dU0). (4.8)I approximate Eqn. (4.7) as mentioned up to linear terms in φ j/U0 and get:dφ jdt= φ j−1(a+c2)−2φ j(a+c2)−φ jv+φ j+1(a+c2)+φ j+1v+c4dd∑k=1(φ j−1+k−φ j−1−k +φ j+1+k−φ j+1−k).(4.9)I can consider stability with respect to each individual Fourier mode, i.e. consider a perturbation ofthe form:φ j = φeλ teθ ji where θ =2pimn, (4.10)n is the total number of lattice sites (or bins) and m = 0, ...,n− 1. For each choice of m there aretwo eigenvalues which leads to a total of 2n eigenvalues. I next substitute this form of perturbationinto Eqn. (4.9) and simplify, using2cosθ = e2pimi/n+ e−2pimi/n,230 1 2 3 4 5 6-6-4-202Figure 4.3: A plot of the real part of the eigenvalues from Eqn (4.11) for three different valuesof c. The onset of instability occurs when Re(λ ) becomes positive (curve crosses thesolid red line). Other parameters used for this plot are: a = 0.5,d = 1 and v = 0.1.to get:λφ = (2a+ c)(cosθ −1)φ + c2d(1+ cosθ − cos(θd)φ − cos((d+1)θ))φ+ v(−1+ cosθ + isinθ)φ .(4.11)These eigenvalues, λ , are complex which is indicative of oscillations. The oscillation period is:2piθ =nm . For perturbations to grow, the real part of the complex eigenvalues needs to be positiveRe(λ ) > 0. Fig. 4.3 shows the graph of Re(λ (θ)) for different values of c. In [9], the results aresimilar but with v= 0. The authors of [9] demonstrated that the onset of instability in their case wasfor θ = 0. In this scenario, the θ = 0 case leads to eigenvalue λ = 0. With this in mind, I assumeinstability first occurs for small value of θ (i.e. choose m to be small).For small θ , the following approximations can be made:sinθ ≈ θ ,cosθ ≈ 1−θ 2/2.(4.12)I use these approximations and simplify the dispersion relation in Eqn. (4.11) to,λ = θ 2(−2a+ cd− v)+ iθv. (4.13)24Thus, (from Re(λ )> 0) I obtain the aggregation condition:cd > (2a+ v). (4.14)As was to be expected, for the case v= 0 the aggregation threshold in Eqn. (4.2) from [9] is retrieved.Notably, this inequality suggests that the addition of a directional bias term makes aggregation lessfeasible. Introducing the directionality, v, may be causing a competition between moving towards atarget to the left (i.e. due to taxis) and moving towards the higher density of cells. As in [9], cd is thenon-dimensional rate after after choosing a direction based on cell density times the sensing range(in number of bins). In order for cells to aggregate, the cd term needs to overcome the persistence inmotion, a, plus the directional bias, v. This result also suggests that larger sensing range, d, makesaggregation more feasible.4.5 Simulation resultsTo compare the instability threshold calculated in Eqn. (4.14) to the discrete model, Eqn. (4.6), Iconduct simulations using Matlab. To avoid boundary effects, all simulations are conducted withperiodic boundary conditions. A small perturbation around the uniform steady state is introducedin the initial condition,u(x,0) = 1+0.01sin(2pix/n)−0.02sin(8pix/n)),where the number of lattice sites (bins) is n = 50 and the size of each one is ∆x = 1. The time stepused is ∆t = 0.1. Choosing a smaller ∆t does not alter the results. However, if ∆t is too large thesimulations break down (because for stability the solution method, ∆x/∆t, needs to work faster thanthe cell movement). The numerical simulations in Fig. 4.4 demonstrate that tuning the directionalityrate, v, can shift the results between the stability and instability regimes. The dimensionless ratesused in Fig. 4.4 are v = 0 (blue solid line), v = 0.01 (red dotted line), v = 0.1 (yellow dashedline) and v = 0.2 (purple dash-dot line). The simulations also show that the size of the directionalbias term changes the shape of the aggregates. For v = 0.1 the peaks are wider and smaller thenfor v = 0.01. The results indicate a good agreement between the calculated aggregation threshold,Eqn. (4.14), and the numerics.4.6 Continuum limit of the discrete modelIn order to study the pattern formation in the instability regime, I first derive the continuum limitof Eqn. (4.6). I set U j(t) = u(x, t) and U j+k(t) = u(x+ kh, t). Here, h = ∆x represents the inter-bindistance in the continuum model and x is the position of the cell where x j = j∆x. I will expandthe equation in a Taylor series. To simplify the calculations, I use the same approach as in [9] and250 10 20 3002468Figure 4.4: Matlab simulation of the discretized model in Eqn. (4.6) with four different valuesfor the directionality rate, v: blue solid line v= 0, red dotted line v= 0.01, yellow dashedline v= 0.1 and purple dashed-dot line v= 0.2. For small speeds the aggregate shapes aresimilar to the results without a directional bias. When the rate v exceeds a set thresholdthe results move from the instability regime into the uniform stable steady-state. Otherparameters in these simulations are: d = 1.0, a = 0.5 and c = 1.2. Results are shown att=10,000 (dimensionless).define:g(h) = u(x−h) ∑dk=1 u(x−h+ kh)∑dk=1(u(x−h− kh)+u(x−h+ kh)), (4.15)to describe the probability of cells choosing a direction based on surrounding cell density. Thenonlinear terms in Eqn. (4.6) then become:U j−1η+j−1+U j+1η−j+1 = g(h)+g(−h)= 2g(0)+h2g′′(0)+h412g(iv)(0)+O(h6).(4.16)Because of the symmetry of the η-terms in the model, the odd terms cancel out. After a longcalculation, carefully keeping track of terms of the same order, I obtaing(0) =u2, g′′(0) =−d2uxx and (4.17a)g(iv)(0) =[−12− d2(d2+2d+3)]uxxxx+12(2d+1)(d+1)2(uxuxxu)x. (4.17b)26To study the stability of the continuum system, it is sufficient for me to keep terms up to orderh2, O(h2). Later, I will investigate the pattern formation of the model in which case I will need toconsider higher order terms because the O(h2) equation becomes ill-posed in the patterning regime.The continuum limit of Eqn. (4.6) up to order h2 (using Eqn. (4.17a)) is obtained by substitutingthe continuous version of the cell density and the calculated g(0) and g′′(0) into Eqn. (4.6):ut = u(x+h, t)(a+ v)+u(x, t)(−2a− c− v)+u(x−h, t)a+ c(g(h)+g(−h))= h2auxx−u(x, t)(c+ v)+ vu(x+h, t)+ c(u(x, t)− h2d2uxx)= vhux(x, t)+h22(2a− cd+ v)uxx(x, t).(4.18)Recall that all the rates in my model depend on the inter-bin distance, h, and thus the h remainsin the PDE here and in the following section. Interestingly, introducing the directional bias in thisasymmetric manner makes it appear in the effective diffusion coefficient.4.6.1 Stability conditionThe PDE in Eqn. (4.18) has a homogeneous steady state: u(x, t) = u0, where u0 is any constant(non-negative for biological relevance) and a non-homogenous steady state. I calculated the non-homogenous steady state by setting the time derivative equal to zero, integrating twice and usingthe boundary condition:ux(x, t) =−h(2a+ v− cd)2v uxx(x, t)⇒ u(x, t) =−h(2a+ v− cd)2vux(x, t)+K1⇒ uss(x, t) = K1e2vx/(h(cd−2a−v))+K2.In the instability parameter regime (from Eqn. (4.14)) the exponent is positive. Thus, the non-homogeneous steady state will blow up for large x (unbounded) which does not make sense for abiological system.One way to compare the continuum limit equation to the discrete model is to do a steady stateanalysis and make sure to retrieve the stability condition from Eqn. (4.14). By introducing thefollowing perturbation around the steady state,u(x, t) = u0+ uˆeλ t+iqx,27I obtain the eigenvaluesλ = ivhq−q2 h22(2a− cd+ v).These are complex eigenvalues andRe(λ )> 0 if cd > (2a+ v).This instability condition is the same as I obtained for the discrete model in Eqn. (4.14), as expected.Note that in the instability regime any wavenumber, q, will lead to instability.4.6.2 Aggregation shapeThe approximation up to O(h2) gave the above instability condition. However, in the regime whereinstability occurs and aggregation is possible, Eqn. (4.18) has anti-diffusion and the solutions blowup. Thus, I need to write the expansion of the continuum system up to O(h4). The expectation is thatthis will result in dispersion term(s) that can dominate the anti-diffusion. I include the calculationfrom Eqn. (4.17b) and obtain:ut = a(h2uxx+h412uxxxx)+ v(hux(x, t)+h22uxx+h36uxxx+h424uxxxx)− h22(cd)uxx(x, t)− ch424(1+d(d2+2d+3))uxxxx+c(2d+1)(d+1)22(uxuxxu)x= vhux(x, t)−Auxx(x, t)+Buxxx−Cuxxxx+D(uxuxxu)x.(4.19)Where the coefficients have been defined as follows:A =h22(cd−2a− v),B =h36v,C =h424(c(1+d(d2+2d+3))−2a− v) andD =h424c(2d+1)(d+1)2.All the parameters A,B,C and D are dimensionless and proportional to the inter-bin size h to somepower (hp). Recall that the dimensionless rates can also depend on h. Therefore, the terms thatsurvive once the limit h→ 0 is taken will vary based on the depends of the parameters on h. (Forexample, if v∼ h−2 then v will remain in A as h→ 0 but not in B and C.)284.6.3 Numerical analysisIn the instability regime and with v > 0, all the aggregate parameters are positive, A,B,C,D > 0. Itcan also be shown that D−C is positive in the instability regime which is important when carryingout a dispersion analysis for the stability of the higher order modes. Interestingly, the effectivediffusion coefficient, A, depends not only on the parameter a, as one might naively expect basedon the cell movement described in the discrete system but also on other quantities, namely: thedirectional bias rate v, the sensing radius d and the rate at which the cell moves after choosing adirection based on cell density c.I used a semi-implicit finite difference scheme to solve Eqn. (4.19) and compare the results tothe discrete model (the same method was used in [9]). Space is discretized with ∆x = h = 1 andu(x, t)∼~u(t) where~u is the vector approximation for u.~u(t+∆t)−~u(t)dt= vhL1~u(t)−AL2~u(t+∆t)+BL3~u(t)−CL4~u(t)+DL1(L1~u(t)L2~u(t)~u(t)).The linear derivative operators, Ln, are approximated using central finite differences:L1u =u j+1−u j−12∆x, L2u =u j+1−2u j +u j−1(∆x)2,L3u =u j+2−2u j+1+2u j−1−u j−22(∆x)3andL4u =u j+2−4u j+1+6u j−4u j−1+u j−2(∆x)4.(4.20)The results are shown in Fig. 4.5. There is good qualitative agreement between the discrete and con-tinuum models. Initially, four aggregates emerge at the same location in both simulations. Aroundt=2,500 the two aggregates on the far left and the two aggregates on the far right merge in the contin-uum model whereas the two aggregates in the middle and the two on the ends merge in the discretemodel. I also conducted simulations with the discrete and continuum models, outside the instabilityregime (same parameters except c= 1.8). In both simulations, the perturbations decay to the steadystate as expected (results not shown). It should be noted that there seems to exist a metastable stateof more than one peak in the system. Similar behaviour has been observed and studied in othermodels such as in Cahn-Hilliard PDEs [79].Next, I rescale the equation using x= xˆ√CA and t = tˆCA2 . xˆ and tˆ are unit-less variables, the factor√CA scales with h andCA2 scales with dimensionless time (rescaled initially by the characteristic timescale). After rescaling (and dropping the hat notation), Eqn. (4.19) becomes:ut = vh√CA3ux(x, t)−uxx(x, t)+√B2ACuxxx−uxxxx+ DC(uxuxxu)x. (4.21)29A) B)0 20 4001000200030001230 20 400100020003000123Figure 4.5: Matlab simulation of the A) discrete model in Eqn. (4.6) and B) continuous modelin Eqn. (4.19). Qualitatively the results are in good agreement. The aggregates mergeat the same time, even though different aggregation peaks merge. The dimensionlessparameters used in both simulations are: a = 1.5, c = 3.6, d = 1 and v = 0.0, n = 50.The initial condition is u(x,0) = 1+0.1sin(8pix/n).A) B)0 20 400100020003000 0.511.50 20 4001000200030000.9511.05Figure 4.6: Matlab simulations using the continuous model, Eqn. (4.19) with non-zero direc-tional bias. A) travelling wave solution for v= 0.05. B) for large bias, v= 0.6, the initialperturbation decays to the uniform steady state. Other parameters and initial conditionsare as in Fig. 4.5.30It can be shown that DC >1 when cd > (2a− v) (the instability condition from Eqn. (4.14)). Forsimplicity, I defineα =DC=c(2d+1)(d+1)2c(1+d(d2+2d+2))+(cd−2a− v) ,β = vh√CA3=v√3(c(1+d(d2+2d+2))+(cd−2a− v)(cd−2a− v)3)1/2,γ =√B2AC=2v√3[(cd−2a− v)((cd−2a− v)+ c(1+d(d2+2d+2)))]1/2,(all unit-less parameters and independent of h). Now, Eqn. (4.21) can be written as:ut = βux(x, t)−uxx(x, t)+ γuxxx−uxxxx+α(uxuxxu)x. (4.22)The addition of the directional bias term results in two extra terms in the PDE as comparedto Eqn. (4.3), namely the first order advection term and the third order term. This PDE is notstraightforward to analyze. It has up to fourth order derivatives and a nonlinear term. However, Iwill study a special solution of the equation.4.6.4 Slow speed regime: asymptotic approximationThe simulations, using both the discrete model in Fig. 4.4 and the continuous model in Fig. 4.6,indicate that the solution looks like a periodic travelling wave. The parameter, v, can be found inall of the coefficients α , β and γ . Therefore, the wave speed is not the advection coefficient β andchanging the reference frame to u(x, t) = u(x− v0t) is not advantageous. This is not trivial, but forthe special case when the speed is slow, asymptotic expansions can be used to find the speed of thetravelling wave in terms of the parameters.To use an asymptotic expansion the parameters need to be chosen such that β and γ in Eqn.(4.22) are small and thus I can writeβ = β0ε, γ = γ0ε, where ε 1.Now the advection and the third order derivative terms are small relative to the rest of the terms inthe equation.I seek solutions of the formu(x, t) = ω (x− εc0t)+ εω1(x− εc0t)+higher powers of εwhere c0 is the wave speed to be determined and ω is the steady state solution of the equation whenε = 0.31For the O(1) problem, I retrieve the steady state of Eqn. (4.3) from [9] where v = 0, namely0 =−ωyy−ωyyyy+α(ωyωyyω)y.Next, from the O(ε) problem, I obtain:− c0ωy = β0ωy+ γ0ωyyy+Lω1, (4.23)where L is the associated linear operator:Lφ :=−φyy−φyyyy+α(φyωyyω+ωyφyyω− ωyωyyω2φ)y.From the Fredholm alternative, I can find a solvability condition for Eqn. (4.23). The Fredholmalternative states that there exists a solution to Lω1 =ω0 if and only if the inner product<ω0,ψ0 >=0 where ψ0 is in the kernel of L∗ (the adjoint operator), that is L∗ψ0 = 0. The adjoint of the operatorL (defined as < ψ,Lφ >=< L∗ψ,φ >) can be found from:∫ψ(φyωyyω+ωyφyyω− ωyωyyω2φ)ydy =∫ (−ψyφyωyyω −ψyωyφyyω+ψyωyωyyω2φ)dy=∫φ((ψyωyyω)y−(ωyψyω)yy+ψyωyωyyω2)dy.(4.24)Thus, I define the operator L∗ asL∗ψ =−ψyy−ψyyyy+α((ψyωyyω)y−(ωyψyω)yy+ψyωyωyyω2).Then, by multiplying Eqn. (4.23) by ψ0 and integrating I obtain:(−c0−β0)∫ψ0ωydy = γ0∫ψ0ωyyydy.This gives the speed in terms of ψ0 :c0 =−γ0∫ψ0ωyyydy∫ψ0ωydy−β0 := γ0I−β0. (4.25)ψ0 is independent of γ0 and β0, it only depends on α. I can now write the speed, c0 in terms of theoriginal parameters and the integral term, I:c0 =−v K(2I+1)+ c(1+G)[3K3(K+ c(1+G))]1/2(4.26)32where K := cd−2a− v > 0 in the instability regime, and G := d(d2 +2d+2)≥ 5 because d ≥ 1.Therefore, the wave speed is negative (c0 < 0) in the instability regime.It is not clear how to obtain ψ0 analytically but it is computable numerically. To do so, Idiscretize L∗ and solve L∗ψ = λψ where ψ0 is the eigenvector corresponding to the eigenvalueλ = 0. Note that Lωy = 0, that means that ψ0 is the adjoint of ωy. However, L does not appear to beself-adjoint (< ψ,Lφ >6=< Lψ,φ >) so in principle ψ0 is not ωy.To obtain the speed numerically, I will start by expanding L∗ :L∗ψ =[3αωωxωxx− 2αω3 (ω3x )]L1ψ+[2αω2ω2x −αωωxx−1]L2ψ−[αωωx]L3ψ−L4ψ:= (L1K1+L2K2+L3K3−L4)ψ.Ln represents the nth derivative operator (as described in Eqn. (4.20)) and Kn are diagonal matriceswithK1( j, j) =3αω( j)ω( j)xωxx( j)− 2αω( j)3 (ωx( j)3)and so forth. From [9] and Eqn. (4.5) I haveω =C[sech√α−12(y+ εc0t)] 2α−1, (4.27)ωy = ω(tanh√α−12(y+ εc0t))( −1√α−1),ωyy = ω(tanh√α−12(y+ εc0t))2(−12−(α−32α−2)).To calculate the wave speed, c0, I solve Eqn. (4.25) using Matlab and I found good agreementwith the simulations of the discrete model, see Fig. 4.7. The simulation results (red dotted line) andthe analytical result (blue solid line) are in good agreement for small directional bias (up to aboutv = 0.4). For higher values of v the analytical solution estimates a higher speed for the aggregatemotion than obtained in the simulations (the estimate gets even worse for v > 0.6, not shown in thefigure). Recall that the asymptotic approximation only holds for small values of β and γ .It is worth noting here that the wave speed in the simulations does not seem to change afterpeaks merge (excluding the short transition period). That suggests that the size of the aggregatedoes not alter the speed at which it moves.Fig. 4.8 contains simulation results of the continuum model, Eqn. (4.19). I used fewer gridpoints, n=30, to get the cells to merge to one aggregate within a reasonable time scale. In thissimulation, the directional bias is small, v=0.01. Fig. 4.9 shows a comparison of the predicted peak330 0.1 0.2 0.3 0.4 0.5 0.600.20.40.60.81Figure 4.7: A comparison of the analytical solution (blue solid line) from Eqn. (4.25) and thesimulation results (red dotted line) of the travelling wave speed for different values of thedirectional bias rate v. Parameters: a = 0.8, c = 3.5, d = 1 and n = 30.shape from [9] (Eqn. (4.27)) to a simulation. The results demonstrate that for this set of parametersthe peak shape derived in the absence of directional bias is a good approximation of the peak shapefor small v.0 10 2005000100001234560 10 20440045004600470048004900 123456Figure 4.8: This simulation is of the continuum model approximation, Eqn. (4.19). A) Kymo-graph of the simulation result with ∆t =0.1. The apparent discontinuity around t = 4700is caused by the two peaks merging rapidly. B) Same kymograph as A) except with thetime interval 4400-5000 and ∆t =0.01, to show the rapid peak merger. The dimensionlessparameters used in this simulation are: a=1.0, c=3.0, d=1.0, v=0.01and ∆x =1.0 and ∆t=0.1. I used n=30 and the initial condition u = 1.0+0.01sin(2pix/n)−0.02sin(8pix/n).340 10 20 300246810Figure 4.9: To the leading order the peak shape takes on the form in Eqn. (4.27) (bluesolid line) which is in good agreement with the simulation result using Eqn. (4.6)(red circles) (the y-axis in the plot represents the dimensionless density of U). Theparameters from Fig. 4.8 are used except a=0.8 and c=3.5. The initial condition isu = 1.0+0.01sin(6pix/n)−0.02sin(2pix/n). The results are shown at t=110,000.4.7 Alternative models and approaches for future work4.7.1 Symmetric model with directional biasAs mentioned previously, the model in Eqn. (4.6) is not symmetric. An alternative, symmetricmodel with directional bias to the left could be derived where the bias increases the probabilityof moving to the left and decrease the probability of moving to the right. In the context of cellmotion, this might for example represent the re-polarization of cells that were moving away from achemotactic source. This model could be of the following form:dU jdt=U j−1(a− v2+ cη+j−1)+U j (−2a− c)+U j+1(a+v2+ cη−j+1). (4.28)35By carrying out equivalent steps as in section 4.6 to obtain a PDE for the model, I obtain the sameform of the PDE as in Eqn. (4.19) but the coefficients are:A =h22(cd−2a) ,B =h36v,C =h424(c(1+d(d2+2d+3))−2a)andD =h424c(2d+1)(d+1)2.The differences between the PDEs for the models in Eqs. (4.6) and (4.28) are in the diffusioncoefficient, A, and the coefficient for the fourth order spatial derivative, C. These coefficients areindependent of the directional bias, v, in the symmetric model. This interesting result would benefitfrom further investigation in future work.4.7.2 Deriving the microscopic parameters from macroscopic observationsIn this chapter, I have derived a discrete model based on movement of individual cells to study how acollection of cells would behave. Alternatively, one can make measurements of how a collection ofcells migrates and try and deduce how the individual cells in the collection behave. If the migrationpattern of the collective can be approximated with a PDE of the form:ut = vux(x, t)+Duxx(x, t)+M1uxxx+M2uxxxx+M3(uxuxxu)x, (4.29)then it can be compared to transition probabilities for individual cells. The transition probabilitiescan be written as:U j(t+∆t) =U j(t)−U j(kR+ kL+ cη+j + cη−j )+U j+1(kL+ cη−j+1)+U j−1(kR+ cη+j−1) (4.30)I use Taylor expansion up to order h2 and obtain:U j(t+∆t)−U j(t) = ∆tU jt =U jxh(kL− kR)+U jxx h22(kL+ kR− cd) (4.31)Now I can compare this to an advection-diffusion equation (first two terms in Eqn. (4.29)) to matchthe unknown parameters dependance on h. I get:kR =D∆th2+cd2− v∆t2h, andkL =D∆th2+cd2+v∆t2h.36Next, I keep higher order terms in the expansion to find c, namely:∆tU jt =U jxh(kL− kR)+U jxx h22(kL+ kR− cd)+U jxxx h36(kL− kR)+U jxxxxh424(kL+ kR− c(1+d(d2+2d+3)))+(U jxU jxxU j)xh424(c(2d+1)(d+1)2).Finally, I set the coefficient of the fourth order derivative equal to M2 and find thatc =11+d(d2+2d+2)(2D∆th2− 24∆tM2h4). (4.32)In this formalism the third order derivative in Eqn. (4.29) disappears when h→ 0 and onlythe first order derivative is a function of v. In future research, it would be interesting to compareexperimental observation of a cell collective to this framework and deduce how individual cellsbehave.4.8 DiscussionIn this chapter, I studied the effects of cell-cell attraction combined with a directional bias term oncell migration and aggregation. I derived a continuum limit of the discrete model and demonstratedthe relationship between the parameters of the models. Interestingly, the diffusion coefficient isa function of all of the parameters of the discrete model which may be counterintuitive (researchquestion 5). This effective diffusion coefficient is reminiscent of the well studied phenomenon ofTaylor dispersion [91]. Often, classical PDE models are formulated with standard phenomenolog-ical terms (i.e. chemotaxis and random motion) to describe behaviour of biological systems. Myresults emphasize that by doing so one needs to be vigilant when interpreting the parameters interms of a biological system. Furthermore, the rates in my model can depend on the size of thebins (h = ∆x) used to write down the transition rates. Thus, quantitative comparison of models andexperiments requires a careful calibration of experimental measurements.The use of a directional bias term in my model can be interpreted as a kind of taxis. One examplewould be chemotaxis where the cells migrate up a gradient of a diffusible chemical. The additionof this directional bias term makes aggregation of cells less feasible, see Eqn. (4.14) (researchquestions 1 and 2). Put differently, the attraction of cells to other cells needs to be stronger when adirectional bias term is added to the equation for aggregation to occur (research question 3). I foundgood agreement between the discrete and continuum models, see Fig. 4.5. However, I was not ableto calculate a general analytic solution for the wave speed of this seemingly simple model. Instead,I calculated the wave speed for the special case when the speed is small and thus the advection termand the third order derivative are small relative to the diffusion and the fourth order term, Eqn. (4.25)37(research question 4). In the small speed limit, the aggregation shape can be approximated as theshape obtained without directional bias from [9], see Fig. 4.9. The speed of the aggregate did notseem to depend on the density of cells, i.e. the speed did not change after two peaks merged.This model can be used to capture an in vitro laboratory experiment with one cell type and somesort of external stimuli. For example, it can model the migration of the cancer cells in Fig. 4.1 thatare attracted to an external source of HGF as well as a chemical signal, CSF-1, from other cancercells. Furthermore, the sensing range can be interpreted as the length over which cells communicateusing protrusions, nanotubes or chemical signals. A typical diffusive length scale (√Dr, where Dis the diffusion coefficient of a chemical signal and r a chemical rate constant) is one example of alength scale based on chemical signalling. Alternatively for contact interactions, the sensing rangeshould be estimated based on the length of protrusions in terms of cell diameters. Naturally the sizeof the sensing range could vary between cell types.The model derived in this chapter is simple and generic, however it can capture basic cell mi-gration properties. The migration rules of the discrete model capture the phenomenological cellbehaviour but physical properties of cells have not been accounted for. Those properties includecell volume, cell shape, membrane tension or any force generation. In other modelling work usingdifferential equations, adhesion or other cell-cell interactions have sometimes been described usingchemotaxis to a chemical signal, integro-PDEs or interaction kernels [41, 68, 69]. Additionally, vol-ume effects have been accounted for by for example allowing cells to only move into empty spaceor by displacing other cells [40].In future work, some advancements can be incorporated into this framework. First, the apparentmetastability of more than one aggregate as well as the stability of a single aggregate should beanalyzed. Second, the previously mentioned volume exclusion effects should be explored. Third,it may be beneficial to consider different analytical ways to study the continuum model. Lastly, arigorous comparison of the asymmetric model in Eqn. (4.6) and the symmetric model in Eqn. (4.28)should be conducted.38Chapter 5Discrete model of interactions betweendistinct cell types: analysis5.1 IntroductionIn this chapter, I will continue to study interactions of cells in the discrete model framework from[9]. I will focus on interactions between two different types of cells. This research is motivatedby the study of chemical signalling between cells that are attracted to a chemical released by cellsof the same type, and cells attracted to a chemical released by a different cell type. One possiblebiological interpretation for such a model would be a short-ranged signalling-loop between breastcancer cells and cells of the immune system, see Fig. 5.1. These interactions enable tumor cells toutilize the immune system to enhance their migration properties. Here, I will not consider the detailsof the chemical signalling but consider effective ranges of interactions (this can be interpreted asapproximating reaction-diffusion equation describing chemical dynamics with a smooth profile).I will derive a simple one-dimensional (1D) model for the cell interactions to explore how theinteractions between cells of the same type and cells of different types alter the aggregation patterns.Using the same analytical framework as in the previous chapter along with numerical simulations, Iwill study the dynamics of the model.5.2 Research questionsI extend the framework from Chavy-Waddy and Kolokolnikov [9], to study interactions betweentwo cell types. The research questions that I will explore in this chapter are:1. Under what parameters do the two cell types form aggregates?2. How do interactions within species (same type of cells) and between species (different typesof cells) affect aggregation?39CSF-1EGFMacro- phageTumor cellBlood vesselHGFEndothelial cellsCSF-1EGFFigure 5.1: Proposed interactions between tumor cells, macrophages (cells of the immune sys-tem) and endothelial cells (cells that line the blood vessels). The endothelial cells releasea chemical signal, HGF, which the tumor cells can detect and migrate towards. Fur-thermore, the tumor cells interact with macrophages via a short-ranged signalling loop(EGF/CSF-1), which enhances the directional migration of tumor cells. Lastly, tumorcells are attracted to one another via a chemical signal (CSF-1/CSF-1R).5.3 Model derivationI consider two different cell types, U and W . As in to the previous chapter, I describe the motilityproperties of the density of these two cell types between bins ( j = 1, ...,n). The dimensionlessrates at which the cells U and W migrate are a1 and a2, respectively. c11 is the migration rate forcells U after choosing a direction towards higher density of cell type U . This can be interpretedas the attraction of cells to other cells of the same type for example through chemical signalling.c12 is the migration rate for cells U after choosing a direction towards higher density of cell typeW . Attraction of cells to another cell type is comparable to paracrine effects (short ranged chemicalsignalling loop between cell types) in biological systems. Similarly, c22 is the migration rate for Wcells going towards higher density of W cells and c21 is the migration rate for W cells switchingdirection towards U cells. The resulting non-dimensional equations for the cell densities in each binare:dU j/dt =U j+1(a1+ c11ηu−j+1+ c12ηw−j+1)+U j(−2a1− c11− c12)+U j−1(a1+ c11ηu+j−1+ c12ηw+j−1),dWj/dt =Wj+1(a2+ c22ηw−j+1+ c21ηu−j+1)+Wj(−2a2− c22− c21)+Wj−1(a2+ c22ηw+j−1+ c21ηu+j−1).(5.1)40The probabilities of cells choosing a directions towards higher density of either U or W cells are:ηu±j =∑dk=1U j±k∑dk=1(U j+k +U j−k), (5.2a)ηw±j =∑dk=1Wj±k∑dk=1(Wj+k +Wj−k), (5.2b)where d is the sensing range of the cells. Here, I have assumed that the sensing range for both celltypes is the same. In general, this range may differ between cell types (see discussion of what candetermine this range in the previous chapter). The η-equations capture the fraction of each cell typewithin d bins of bin j that are to the left (-) or right (+).This system has a constant uniform steady state, U0 ≥ 0 and W0 ≥ 0. I am interested in derivinga condition for the appearance of a non-homogeneous steady state.5.4 Steady-state analysisI use linear stability analysis to study the stability of the steady-state of this model. I introduce smallperturbations, |φ j| U0 and |ψ j| W0, around the steady states, namelyU j =U0+φ j and Wj =W0+ψ j.Eqn. (5.1) for U can then be written in the following way:dφ jdt= (U0+φ j+1)(a1+ c11ηu−j+1+ c12ηw−j+1)+(U0+φ j)(−2a1− c11− c12)+(U0+φ j−1)(a1+ c11ηu+j−1+ c12ηw+j−1),(5.3)where the η-terms are of the form:ηu+j−1 =∑dk=1(U0+φ j−1+k)∑dk=1(U0+φ j−1+k +U0+φ j−1−k)). (5.4)By using Taylor expansion and only keeping terms that are linear in φ j/U0, I obtain:ηu+j−1 ≈12+14dU0(d∑k=1(φ j−1+k−φ j−1−k)). (5.5)41I insert this linear approximation into Eqn. (5.3) and get:dφ jdt= (φ j−1+φ j+1)(a1+ c11/2+ c12/2)−φ j(2a1+ c11+ c12)+c114d(φ j+d+1+φ j+d−φ j−φ j+1+φ j−d−1−φ j−d−φ j−φ j−1)+c124dU0W0(ψ j+d+1+ψ j+d−ψ j−ψ j+1+ψ j−d−1−ψ j−d−ψ j−ψ j−1).(5.6)This is a (2n)× (2n) linear problem which can be decomposed into n 2×2 problems by using thesame type of perturbation as in Eqn. (4.10):φ j = φeλ teθ ji and ψ j = ψeλ teθ ji where θ =2pimn,where n is the total number of bins and m = 0, ...,n− 1. Therefore, Eqn. (5.6) can now be writtenas:λφ = φ(cosθ −1)(2a1+ c11+ c12)+(c114dφ +c124dU0W0ψ)(2cos((d+1)θ)+2cos(dθ)−2−2cosθ).(5.7)For small θ (nm), I use the approximation in Eqn. (4.12), in which case, Eqn. (5.7) simplifies to:λφ = φ−θ 22(2a1+ c11+ c12)+(c112dφ +c122dU0W0ψ)dθ 2(1+d). (5.8)Following the same procedure, now for W gives,λψ = ψ−θ 22(2a2+ c22+ c21)+(c222dφ +c212dW0U0ψ)dθ 2(1+d). (5.9)The Jacobian for this set of equations is:M =θ 22[−(2a1− c11d+ c12) c122 U0W0 (1+d)c212W0U0(1+d) −(2a2− c22d+ c21)]. (5.10)For simplicity, I will first consider a special case, i.e. c11 = c21 = 0 and d = 1. In this scenario, bothcell types are only attracted to the cell type W . Then:M∗ =θ 22[−(2a1+ c12) c12U0W00 −(2a2+ c22)+2c22]. (5.11)42The trace (Tr) and determinant (Det) of this matrix are,Tr(M∗) =−2(a1+a2)− c12+ c22 andDet(M∗) = (2a1+ c12)(2a2− c22).In this example, I can never have Det(M∗)>0 and Tr(M∗)> 0. Thus, instability only emerges whenDet(M∗)< 0, which occurs for:c22 > 2a2.This instability threshold matches simulation results, see Fig. 5.2. Intuitively, this threshold seemsreasonable since the cell type U does not attract other cells and thus I expect this threshold to beequivalent to the one obtained in [9] (Eqn. (4.2)), as is the case.Now, getting back to the full system, I have the following trace and determinant:Tr(M) =−(2a1+ c12− c11d)− (2a2+ c21− c22d) andDet(M) = (2a1+ c12− c11d)(2a2+ c21− c22d)− c12c214 (1+d)2.The discriminant for this system is:Tr(M)2−4Det(M) =((2a1+ c12− c11d)− (2a2+ c21− c22d))2+2c11c22d2+ c12c21(1+d)2 ≥ 0.Since I have assumed that all the ci j parameters are positive, the discriminant is always positive andthe system only has real eigenvalues, no Hopf bifurcation or oscillations. (However, for the specialcase where c12 and c21 have opposite sign, there is a Hopf bifurcation and oscillations occur in thedensity of the two species. I will explore this numerically in the coming section.)The onset of instability requires (i) Det(M)< 0 and thus,c12c214(1+d)2 > (2a1+ c12− c11d)(2a2+ c21− c22d), (5.12)alternatively, (ii) Det(M)>0 and Tr(M)> 0 and therefore instability occurs when,c22d > 2a2+ c21, c11d > 2a1+ c12 and(2a1+ c12− c11d)(2a2+ c21− c22d)> c12c214 (1+d)2.(5.13)The parameters c11 and c22 represent attraction to the same type of cell (homotypic) and c12 andc21 represent attraction to a different type of cell (heterotypic). In (i), it suffices to have the homo-typic attraction of one of the cells strong enough to overcome the persistence in motion ai and theheterotypic attraction. Alternatively, if both heterotypic attractions are none zero, their strengths as43well as a large sensing range, d, can cause instability. In (ii), the homotypic attraction for both celltypes is strong again (large c11 and c22). In this case, I have an additional condition which puts anadded constraint on the magnitude of the heterotypic terms, c21 and c12. However, instability stilloccurs if this constraint is violated because of the condition in (i) (so long as Det(M) 6=0).5.5 Continuum limit of the discrete modelI next derive a continuum limit for this model to explore aggregation patterns. I set U j(t) = u(x, t),U j+k(t) = u(x+ kh, t), Wj(t) = w(x, t) and Wj+k(t) = w(x+ kh, t) where h = ∆x is the inter bindistance. Using functions g1 and g2 for the non-linear terms as in the previous chapter (Eqn. (4.15)),I define:g1(h) = u(x−h) ∑dk=1 u(x−h+ kh)∑dk=1(u(x−h+ kh)+u(x−h+ kh))andg2(h) = u(x−h) ∑dk=1 w(x−h+ kh)∑dk=1(w(x−h+ kh)+w(x−h+ kh)).(5.14)Notice that g2 is a function of both the cell types u and w because it captures the attraction of celltype u to higher densities of cell type w. The nonlinear terms for U in Eqn. (5.1) then become:U j−1η1+j−1+U j+1η1−j+1 = g1(h)+g1(−h) = 2g1(0)+h2g′′1(0)+h412g(iv)1 (0)+O(h6),U j−1η2+j−1+U j+1η2−j+1 = g2(h)+g2(−h) = 2g2(0)+h2g′′2(0)+h412g(iv)2 (0)+O(h6).(5.15)As in the previous chapter, I haveg1(0) =u2and g′′1(0) =−duxx2.By carrying out equivalant calculations for g2(h), I getg2(0) =u2and g′′2(0) =uxx2− uwxxwd+12+w2xuw2(d+1)2− uxwxwd+12.Since g2 is a function of both u and w, I no longer have as many terms cancelling in the calculationsfor g′′2 (the same is true for higher order derivatives, see Supporting Information A.2). I investigatethe stability of the system by carrying out the approximation up to O(h2). In order to study patternformation, I will need to include higher order terms. The discrete model for U can be approximatedwith the following PDE:44ut = uxxh22(2a1+ c12−dc11)−(h2c12d+12)(wxxuw− (wx)2 uw2 +wxuxw).(5.16)Conducting similar calculations for the cell type w gives the following system of equations:ut = α1uxx−β1(wxxuw− (wx)2 uw2 +wxuxw),wt = α2wxx−β2(uxxwu− (ux)2 wu2 +uxwxu),(5.17)whereα1 =h22(2a1+ c12− c11d), β1 = h2c12 d+12 ,α2 =h22(2a2+ c21− c22d) and β2 = h2c21 d+12 .Recall that the rates in my model are dependent on the bin size h (see comments on this dependanceon h in the previous chapter and Supporting Information A.2 for rescaled version of the equationwith all parameters independent of h).I can perform the following two checks to ensure I can reduce this result to those of the modelfrom [9]. For the bacteria model c12 = c21 = 0, therefore β1 = β2 = 0 and I simply have the diffusionequation. Or, I can set u=w in which case the first order derivative terms will cancel and I am againleft with the diffusion equation.5.5.1 Linear stability analysisI can use the above continuum equations to find the instability condition. I introduce the followingperturbations around the uniform steady statesu = u0+ uˆeλ t+iqx and w = w0+ wˆeλ t+iqx.I insert these into Eqn. (5.17) and drop the hat notation for simplicity. I only keep the linear termsand get:uλ =−q2uα1+β1(q2u02w0w),wλ =−q2wα2+β2(q2w02u0u).(5.18)The Jacobian matrix for this system isM =−q2[α1 − u02w0β1− w02u0β2 α2]. (5.19)45The trace and determinant of the matrix M areTr(M) = α1+α2 and Det(M) = (α1α2)− β1β24 .Now substituting the original parameters in for the α’s and β ’s I get:Det(M) =h44[(2a1+ c12− c11d)(2a2+ c21− c22d)− 14(d+1)2c12c21]andTr(M) =−h22(2a1+ c12− c11d+2a2+ c21− c22d) .This is the same instability conditions as I obtained in Eqs. (5.12) and (5.13) for the discrete model,suggesting that this continuum limit is a good approximation of the discrete model.To study pattern formation for this system of equations, I would need to include the next or-der in the approximation because in the instability regime there is anti-diffusion in one or both ofEqs. (5.17). Including those terms turned out to be long and not too informative, see SupportingInformation A.2 for calculations. Instead, I will proceed by solving the discrete system numerically.5.6 Simulation resultsThis modelling framework only allows a limited amount of analytical insights into this system oftwo interacting cell types. Thus, I proceed by studying the system using numerical methods. Iuse periodic boundary conditions, as in the single cell model in the previous chapter. The initialconditions are:U = 1+0.01sin(2pix/n) and W = 1+0.02sin(2pix/n).The number of grid points is n = 50 with spatial resolution ∆x = 1 and time step ∆t = 0.1. As inchapter 4, there seems to exist a metastability of more than one aggregate in this system.5.6.1 The sensing radius and attraction between-species influence the dynamicsIn Fig. 5.2, I conduct simulations with various values for the sensing range, d, and the attraction ofcell type W to cell type U , c21. Increasing these parameters can move the results between the stabil-ity and instability regimes, as is predicted in Eqn. (5.12). Increasing d results in larger aggregatesand the peaks get narrower as can be seen in column two. Increasing the parameter c21 acceleratesthe merging of aggregates into a single aggregate and it alters the ratio between the peaks of the twocell-types. However, increasing the sensing range, d, can make it harder for the peaks to aggregate,third column. Recall that the sensing range does not only affect the attraction to the other cell typebut also to its own type. Therefore, this might be due to the cohesion of aggregates increasing.460 500510d = 1, c21 = 00 500510d = 1, c21 = 0.50 500510d = 1, c21 = 10 500510d = 1, c21 = 20 500510d = 2, c21 = 00 500510d = 2, c21 = 0.50 500510d = 2, c21 = 10 500510d = 2, c21 = 20 500510d = 3, c21 = 00 500510d = 3, c21 = 0.50 500510d = 3, c21 = 10 500510d = 3, c21 = 20 500510d = 4, c21 = 00 500510d = 4, c21 = 0.50 500510d = 4, c21 = 10 500510d = 4, c21 = 2Figure 5.2: Simulation of Eqn. (5.1) for various values of d and c21. Increasing the sensingradius d makes the peaks higher and their profile narrower. Higher values of c21 seem tofacilitate the peaks merging to one. The y-axis represents the dimensionless densities ofU and W . Other parameters in these simulations are a1 = 0.5, a2 = 2.0, c11 = c22 = 0.7and c12 = 0.5. Results are shown at T= 200,000.475.6.2 Including a directional bias for one cell typeI can investigate, numerically, how a directional bias affects the results. In Fig. 5.3, I conductsimulations with bias, v, to the left for the cell type W (equivalent to the implementation in chapter4).dWj/dt =Wj+1(a2+ v+ c22ηv−j+1+ c21ηu−j+1)+Wj(−2a2− v− c22− c21)+Wj−1(a2+ c22ηv+j−1+ c21ηu+j−1).(5.20)Here, the drift increases the probability of moving to the right resulting in an asymmetric equation.For d = 1, increasing the bias to v = 2 moves the results back into the stability regime. This is whatI expected from the results for the model in the previous chapter (Eqn. 4.14). If multiple aggregatesform, the larger d makes it harder for them to merge. The magnitude of the drift, v, changes theshape of the aggregation peaks. Large v, makes it difficult for the cell type U to keep up with thecell type W and thus the aggregate peaks for U decrease as v increases.5.6.3 Biological system: cancer cell-macrophage interactionsAs an example of how this model can be used to study biological systems, I conducted simula-tions that can be interpreted in terms of the tumor cell and macrophage interactions [35, 100], seeFig. 5.4. The tumor cells, cell type U in the model, are attracted to a chemical (CSF-1) released bymacrophages, W, (interpreted as c12 in the model) and the macrophages are attracted to a chemical(EGF) released by tumor cells (c21). Additionally, tumor cells can also be attracted to other tumorcells (c11). Increasing c11 when c12 > 0 increases the peak size for U. However, increasing c12 in thetwo bottom rows decreases the peak size for U. This is to say that the relative strength of attraction tocells of the same type and to cells of different type alters the ratio between the peaks of the two celltypes. This result is in good qualitative agreement with experimental results from [75] where theauthors compared the effect between cancer cells that were only attracted to macrophages, throughparacrine signalling, to those that were attracted to macrophages and tumor cells. The experimentsshow that the ratio between tumor cells and macrophages that migrate together into a micro-needleincreases from 3:1 to 12:1 for cancer cells that are attracted to other cancer cells.5.6.4 Hopf bifurcationAs mentioned earlier, the system can have a Hopf bifurcation, suggestive of oscillatory patterns,when c12 and c21 have opposite signs (meaning that cell type a repels cell type b and cell type b at-tracts cell type a). Fig. 5.5 contains kymographs of simulation results in this parameter regime withc12 = −1.0, cell type W repels cell type U and for four different values of c21. As c21 is increasedthe patterns change from a steady state to travelling waves, with peaks merging and splitting, topulses.480 500510d = 1, v = 00 500510d = 1, v = 0.40 500510d = 1, v = 0.80 500510d = 1, v = 20 500510d = 2, v = 00 500510d = 2, v = 0.40 500510d = 2, v = 0.80 500510d = 2, v = 20 500510d = 3, v = 00 500510d = 3, v = 0.40 500510d = 3, v = 0.80 500510d = 3, v = 20 500510d = 4, v = 00 500510d = 4, v = 0.40 500510d = 4, v = 0.80 500510d = 4, v = 2Figure 5.3: Simulation results of the discrete model with drift, v, to the left for cell type W , asoutlined in Eqn. (5.20) and various values of d. The peaks of cell type W are followed bythe cell type U . The results show how adding drift can move the system into the stabilityregime (top right corner). The y-axis represents the dimensionless densities of U and W .Other parameters in these simulations are a1 = a2 = 0.5, c11 = c21 = 0.0, c12 = 0.4 andc22 = 2.0. Results are shown at t=300,000.490 10 20 300510c11 = 1, c12 = 10 10 20 300510c11 = 1, c12 = 20 10 20 300510c11 = 1, c12 = 40 10 20 300510c11 = 2, c12 = 10 10 20 300510c11 = 2, c12 = 20 10 20 300510c11 = 2, c12 = 40 10 20 300510c11 = 3, c12 = 10 10 20 300510c11 = 3, c12 = 20 10 20 300510c11 = 3, c12 = 4Figure 5.4: Simulations to investigate the affect of both attraction to another cell type (betweenspecies) and attraction to same types of cells (within species). The simulations in thisfigure are conducted with different c11 and c12 parameters. In these simulations, c22 =0, meaning that the cell type W is not attracted to itself. Increasing c11 increases thepeak size for U while increasing c12 decreases the peak size. The y-axis represents thedimensionless densities of U and W . Other parameters in these simulations are a1 = a2 =0.5, c21 = 0.1 and d = 2. Results are shown at t=200,000.50U, c21 = 00 20 4002004006008001000246810W, c21 = 00 20 40020040060080010000.511.522.5U, c21 = 0.50 20 400200400600800100024681012W, c21 = 0.50 20 40020040060080010001234U, c21 = 10 20 4002004006008001000123W, c21 = 10 20 40020040060080010000.511.522.53U, c21 = 20 20 40020040060080010000.511.52W, c21 = 20 20 40020040060080010000.511.52Figure 5.5: Kymographs showing simulations results of the discrete model for the special casewhere one cell type can repel the other. In the first column, c21=0.0 and the systemreaches a steady state. In the remaining columns for c21 > 0 , there are oscillationsin the peaks for the cell density. U : top row. W :bottom row. Other parameters area1 = a2 = 1.0, c11 = 3.0, c22 = 2.5, d = 1 and c12 =−1.0.Chase-and-run In a biological system, this behaviour can occur, for example due to a combinationof chemotaxis of one cell type and chemorepulsion of another or the interplay between adhesion andanti-adhesion. As a particular example, there is the so-called chase-and-run phenomena betweenneural crest cells and placodal cells (these are embryonic cells that interact for example in thedevelopment of the central nervous system). The neural crest cells chemotax towards the placodalcells but the placodal cells run when they are contacted by neural crest cells [92]. This type ofbehaviour can be obtained in a particular parameter regime in our system, see Fig. 5.6. If therepulsion is big, first column, than the peaks split in both directions because of the periodic boundarycondition (the repulsion is felt from both sides of an aggregate). However, for small repulsionvalues, last column, the species W (blue dotted line) follows species U . When W gets close to U itslows down while U speeds up so the distance between the peaks varies. Therefore, W (the placodalcells) chases U (the neural crest cells) and U runs away from W .510 500510T = 2200, c12 = -10 500510T = 2200, c12 = -0.60 500510T = 2200, c12 = -0.30 500510T = 2200, c12 = -0.150 500510T = 2500, c12 = -10 500510T = 2500, c12 = -0.60 500510T = 2500, c12 = -0.30 500510T = 2500, c12 = -0.150 500510T = 2800, c12 = -10 500510T = 2800, c12 = -0.60 500510T = 2800, c12 = -0.30 500510T = 2800, c12 = -0.150 500510T = 3100, c12 = -10 500510T = 3100, c12 = -0.60 500510T = 3100, c12 = -0.30 500510T = 3100, c12 = -0.15Figure 5.6: Simulations of the discrete model showing a parameter regime where a chase-and-run phenomenon is observed. U attracts W whereas W repels U . U : red solid line.W :blue dotted line. The y-axis represents the dimensionless densities of U and W . Otherparameters are a1 = a2 = 1.0, c11 = c22 = 2.5, d = 1 and c21 = 0.1.525.7 DiscussionIn this chapter, I studied interaction between two different cell types that can lead to aggregation.Similarly to the previous chapter, I used a framework to describe these interactions in a discretemodel and derived a continuum limit for the purpose of analysis. I was able to calculate the thresh-old for pattern formation, see the instability condition for the homotypic and heterotypic attractionin Eqs. (5.12) and (5.13) (research questions 1). As one would expect, each of the cell types canaggregate on its own and that instability condition, Eqn. (5.13), is equivalent to the previous re-search with a single species [9]. For aggregates containing both cell types to emerge the heterotypicattraction needs to satisfy Eqn. (5.12) (research question 2). I found that there was good qualitativeagreement between the continuum and discrete models. However, the derivation of a continuummodel to study the patterns that emerge resulted in a PDE that was long and not easily interpreted.Instead, I explored the discrete model using simulations.One example of a biological system where this model can be applied is for interactions betweentumor cells and cells of the immune system called macrophages [35, 75], see Fig. 5.1. Experimentsshowed that when these two cell types were collected from a mouse tumor (using a micro-needle)either with the within species attraction of tumor cells or without it, the ratio between the cell typeswent from 12:1 (tumor cells to macrophages) to 3:1 [75, 100]. These results are in good qualitativeagreement with simulation results in Fig. 5.4. In the top row of the figure, the within species effectfor cell type U was turned off. In the following rows the magnitude of the effect was increased,resulting in the tumor cell-macrophage ratio between the peaks increasing.Another biological example is the chase-and-run phenomena observed between neural crestcells and placodal cells. The results for the interactions between the two cell types having oppositesigns can be seen in Fig. 5.6. The last column of the figure shows how the cell type W chasesthe cell type U . When it catches up the two cell types temporarily stall and then the cell type Uruns away from the cell type W . Therefore, the distance between the peaks of the two cell typesoscillates. If the magnitude of the repulsion from cell type U on cell type W was increased (the firstthree columns in Fig. 5.6) the shape and number of peaks of W varied over time.In future work, the advancements for this framework outlined in the previous chapter (such asincluding physical cell properties), should be investigated. Further analysis of the oscillatory pat-terns that emerge in the parameter regime where a Hopf bifurcation occurs may yield interestinginsights. The model would be most useful if designed to study a specific biological system and totest different hypotheses regarding cell migration regulation. However, I added seemingly smalllevel of details into the model, a second cell type, and found that the analysis in this frameworkquickly got complicated or even impossible. I was able to calculate the regime in which patterningoccurred but beyond that the analysis provided limited insights into migration behaviour. In the fol-lowing chapters, I will focus on how cell-cell interactions influence the dynamics of cell migration.Consequently, I will rely heavily on simulations.53Chapter 6Collective migration in zebrafishdevelopment6.1 IntroductionCell migration plays a key role in many developmental processes. In this chapter, I focus on aspecific developmental aspect of zebrafish development, the migration of the posterior lateral lineprimordium (PLLP). The migrating PLLP deposits neuromasts in its wake that become sensoryhair cells enabling the fish to respond to fluid flow. The zebrafish embryo has the advantages of (a)embodying vertebrate biology (placing it closer to mammalian systems than other model organismssuch as Drosophila or C. elegans) (b) being optically transparent (c) providing a number of cellbehaviours and chemical signalling that can be related to other biological systems, such as cancer.The migrating PLLP is a dynamic structure, with changes in behaviour and cell morphology. ThePLLP consists of more than 100 cells that migrate together, along the length of the zebrafish embryofrom head to tail (approximately 20 hours post-fertilization). The PLLP completes its migration overthe course of a day, migrating around 3 mm [34].The changes in cell behaviour across the PLLP are thought to be determined by a number ofsignalling networks. The PLLP is dominated by Wnt and FGF signalling. Cells at the front (rear)of the PLLP are dominated by WntR (FGFR) activity. [49, 63]. Through a mutually inhibitorysignalling network, Wnt and FGF signalling are responsible for dividing the PLLP into the leadingand trailing zones. Moreover, the Wnt and FGF signalling networks modulate chemokine receptorexpression in the leading and trailing zones, with the receptor CXCR4b expressed largely in theleading zone, and CXCR7b expressed in the trailing zone [2]. The PLLP migrates on a surfaceof cells that express a stripe of the chemokine CXCL12a (also known as SDF-1a) [17, 38, 95]. Inaddition, the leading and trailing cells cooperate to organize migration of the PLLP. The leadingcells are responsible for following chemokine cues, while the trailing cells undergo chemotaxis54toward the FGF ligand produced by the cells in the leading zone [16].In order for a normal developmental sequence to take place, the cell mass in the PLLP must po-larize and select a direction for migration. Failure of any part of the underlying chemical signalling,or disruptions in the mechanical properties, give rise to mutants with specific migration or otherdefects, providing an opportunity to tease apart some of the underlying interactions.Understanding the PLLP from experiments and verbal arguments alone is limited and unlikelyto yield in-depth knowledge. Theoretical analysis and computational cell-based simulations allowthe study of how a given hypothesis or set of hypotheses fare based on experimental observations inthe literature.Section 6.2 will outline the questions that motivate this research. In section 6.3, I will describethe two modelling frameworks that I use in this research, namely a 1D model for the tissue segre-gation as well as 3D simulations of the PLLP migration using individual cells in the HyDiCell3Dframework. Section 6.4 will depict how the theoretical model can lead to the subdivision of cellswithin the cluster and will also include 1D simulations of the theoretical model. Section 6.5 willdemonstrate simulations using the HyDiCell3D to investigate PLLP migration as well as the forcesrequired to maintain cluster cohesion. In section 6.6, I will compare the modelling results to ex-perimental data. Finally, section 6.7 will outline the research conclusions as well as discuss andcompare the results to previously published models.6.2 Research questionsFollowing are the key questions that guide this research:1. Can Wnt-FGF signalling alone give rise to the initial polarization of the primordium intodistinct leading and trailing domains?2. How do the receptor-ligand signalling parameters and cells’ sensitivities to signalling levelsaffect the proportion of the tissue that develops into Wnt and FGF signalling domains?3. Are the proposed interactions with the chemokine CXCL12a sufficient to explain the migra-tion of the primordium?Some of these questions have been partially addressed in previous modelling papers [1, 12, 16,20, 88]. These models have examined the interactions of the PLLP with CXCL12a in migration [88],the migration and neuromast deposition [12, 20], the mechanical conditions for efficient migration[1] and the outcome of laser ablation of parts of the PLLP [16]. All efforts contribute to piecingtogether aspects of the system, though a comprehensive simulation of the entire process is yet tobe carried out. In the Discussion, these previous modelling efforts are compared to this researchproject.556.3 Modelling frameworks6.3.1 Chemical signallingFirst, a simple chemical signalling model is proposed and used to address the initial stage of po-larization of the PLLP into a WntR active front and an FGFR active rear. While the biologicalliterature provides evidence that mutual inhibition between WntR and FGFR activities is involved,it remains to be determined how that inhibition is mediated, and what details control the boundarybetween the two domains.The following assumptions are used in the model design. (1) A given cell is considered to be aWntR (respectively FGFR) active cell if the level of cell surface Wnt (respectively FGF) receptorsis high (relative to some baseline level). (It is possible for both of these to be above a set thresholdin which case the cell is considered WntR). (2) It is known that there are mutually inhibitory inter-actions between FGFR and WntR activities [13]. Here, a receptor-based mechanism is proposed forthis inhibition. Consequently, the assumption is made that bound Wnt (respectively FGF) receptorson the cell surface activate signalling pathways inside the cell that lead to inhibition of expressionof receptors of the other signalling type. Activation of Wnt signalling determines the expression offactors that directly (sef ) or indirectly (dusp6) interfere with efficacy of signaling by FGF receptors.In general, mutual antagonism could be based on mutual interference by one type of bound receptorwith the surface presentation of the other type, see e.g., [2]. As a result, bound FGF receptors on acell contribute to inhibition of Wnt receptors on its surface and vice versa, see Fig. 6.1.6.3.2 Cell-surface receptor levelsThe state of each cell is represented by the number and type of receptors (Wnt vs. FGF) on thecell surface. Cells that have mostly Wnt receptors on their surface are denoted “Wnt expressingcells” and similarly for FGF. In order to allow for cell commitments to evolve with time, mutuallyinhibitory interactions are implemented between Wnt and FGF signalling as described below.Let R(t) be the level of expression of a given receptor type on the surface of a given cell (forexample, in units of nM). Then the rate of change of R(t) is modelled bydRdt= presentation rate− endocytosis rate. (6.1)Both the presentation and endocytosis rates are regulated by signalling networks. Here, the simplestmodel is assumed, in which mutual FGF-Wnt inhibition affects the presentation rate of the antag-onist receptors (as in, for example, [82]), but results are very similar for regulation via dampingof presentation or accentuation of endocytosis rate or both simultaneously (see, e.g. [44], whereanalogous small positive and negative feedback circuits are analyzed in detail). It is also possible toinclude positive feedback of signalling on the rate of presentation of its own receptor type. Except56Figure 6.1: The geometry and signaling of the posterior lateral line primordium (PLLP). Left:A sketch of the PLLP showing side and top-down views, with leading (red) and trailing(green) cells on a stripe of CXCL12a. Right: A schematic diagram of signalling in thePLLP, showing the mutual inhibition between FGFR and WntR signalling. (Black textdenotes model components, white text aides in interpretation of experiments.) FGF sig-naling inhibits Wnt signaling by determining expression of Dkk1; Wnt signaling cellsexpress the gene sef that inhibits FGF signaling. WntR active cells are sources of bothFGF and Wnt ligands. Experimentally, pea3 and lef1 expression levels are used to iden-tify FGF and Wnt signalling, respectively. The WntR-FGFR activity polarization setsup chemokine polarization (CXCR4b vs. CXCR7b). In the model, this leads to the cre-ation of a gradient of CXCL12a that enables directed migration of the PLLP. Left figurecourtesy of Damian Dalle Nogare.for a few specific scenarios (explained in Supporting Information of [47]) the results are largelysimilar.FR(x, t),WR(x, t) are defined to be the concentration of cell surface receptors (in units of nM),and FB(x, t),WB(x, t) to be the concentration of bound receptors at position x. (In general, x denotesa position within the PLLP in the full 3D model. In the 1D model reduction, the average across thewidth and thickness of the PLLP is taken, and explorations are restricted to variations of signallinglevels across the length. In that case, x represents position along the primordium 0≤ x ≤ L, wherex = L is the front of the primordium and x = 0 is the back). The respective FGF and Wnt ligandsare denoted F(x, t),W (x, t).Receptors can be in either unbound (free) or ligand-bound (B) form on the cell surface. Thetotal level of receptors, FR,WR satisfyFR(x, t) = FB(x, t)+Ff ree(x, t), WR(x, t) =WB(x, t)+Wf ree(x, t). (6.2)The mutually inhibitory signalling within a cell is assumed to depend on FB,WB based on the ideathat only the bound receptors signal to downstream intracellular regulatory networks governing57receptor synthesis and presentation.Since ligand binding to receptors is fast on the timescale of the receptor presentation and endo-cytosis rates, the bound receptors are assumed to be in quasi steady state (QSS) with the availableligand. (See Fig. 6.2.) Then, by standard Michaelis-Menten kinetics,FB(x, t)≈ FR(x, t)F(x, t)KF +F(x, t) , WB(x, t)≈WR(x, t)W (x, t)KW +W (x, t), (6.3)where KF ,KW are equilibrium binding constants that are related to the rates of association anddissociation of the ligand-receptor complex (Fig. 6.2 and Eqn. (A.12e) in Supporting InformationA.3.1).*Free ligandFree receptorphosphorylated receptorgene regulationproteinsynthesisendocytocispresentationkonkok2FASTSLOWNUCLEUSFigure 6.2: Ligand-receptor dynamics forms a basis for the model. FGF or Wnt receptorsare synthesized, presented on the cell surface, and, after binding to ligand, phosphory-lated and internalized by endocytosis. Here,the ligand-receptor binding is assumed to berapid on the timescale of such receptor dynamics, so that the bound receptor level can beapproximated by Eqs. (6.3). Figure courtesy of Leah Edelstein-Keshet.Based on Eqn. (6.1) and assumptions above, the evolving levels of Wnt and FGF cell surfacereceptors are modeled as:dFRdt= EF(WB)− rFFR (6.4a)dWRdt= EW (FB)− rWWR (6.4b)where receptor endocytosis rates, rF ,rW are assumed constant, and rates of presentation, EF(WB),EW (FB)of the receptors are modified by mutually inhibitory feedback. Note that the delay between bindingof receptors and inhibition of receptor presentation pathways is neglected. The rates EF(WB),EW (FB)58are decreasing functions of bound antagonist receptors. In particular, simple, commonly used de-creasing Hill functions typical for such inhibitory terms, e.g., [32], are used:EF(WB) =EF,11+(WB/W0)n, EW (FB) =EW,11+(FB/F0)m. (6.4c)Here Ei,1 is the basal rate of receptor presentation when there is no inhibition, and F0,W0 are the“IC50” inhibition parameters with units of bound receptor density. For example, when bound FGFreceptors are at the level FB(x, t) = F0, then the Wnt receptor presentation rate is reduced to 50% ofits basal level. n,m≥ 2 are Hill coefficients. The larger these values, the sharper the switch turningoff the presentation of receptors. Typically n = m = 5 in the model, which slightly increases theregime of bistability of the system compared to the case n,m ≈ 2, but similar behavior is obtainedfor n 6=m and for smaller values of these Hill coefficients. The limit n→∞ leads to analytic insights,as noted in [73]. In brief, negative feedback is proportional to decreasing Hill functions of the formg(c) =11+ cm, 0≤ g(c)≤ 1where c is a dimensionless quantity (ratio of bound receptor level to IC50 parameter). As expected(see e.g. [82]), this type of mutually-inhibitory mechanism can, under appropriate conditions, giverise to subdivision of the tissue into two domains.This Wnt-FGF receptor model is a minimal model that can set up tissue polarization. In othercell types, similar antagonistic effect between Wnt and FGF has been observed [93]. Other modelvariants could include FGF-bound receptors up-regulating Wnt receptor endocytosis, as well asadditional positive feedback (e.g. from FGF bound receptors to FGF receptor presentation) [2] orother signalling cascades. In practice, given the lack of detailed data for such terms, the minimalmodel variant was chosen for the time being. Furthermore, this research focuses on qualitative,rather than quantitative aspects of the model, since parameter values and biochemical details are notavailable.6.3.3 FGF and Wnt ligand concentrationThe ligands diffuse (with diffusion coefficients DF ,DW ), bind to free receptors (at rates kon,i), andhave some turnover rates δi. When the ligand concentrations are normalized by baseline (steadystate) values, the dimensionless equations for ligand concentrations are:∂F∂ t=DF∂ 2F∂x2−δFF−κFFR F1+F +ρFWRW1+W(6.5a)∂W∂ t=DW∂ 2W∂x2−δWW −κWWR W1+W +ρWWRW1+W(6.5b)59whereDi =DiL2, κF = kF,onF1, κW = kW,onW1, (6.5c)andρF =pFW1KF, ρW =pWW1KW. (6.5d)(See details of scaling in Supporting Information A.3.1.) Here Di =DF ,DW are rates of diffusion ofthe FGF and Wnt ligands, δF ,δW are the rates of decay of the ligands due to nonspecific degradationin the PLLP, ki,on are rates of removal by binding to ligand-specific receptors of the given type andF1,W1 are steady state receptor levels when there is no negative feedback. The terms, ρFWR W1+Wand ρWWR W1+W , in Eqs. (6.5) model the production of FGF and Wnt ligands respectively by WntRactive cells. (See the derivation of these terms in Supporting Information A.3.1.)It is to be emphasized, that the model is a qualitative, rather than quantitative predictive modelsince much of the biological detail, and most of the biological rate constants are unknown and/or noteasily measurable. Nevertheless, an attempt has been made to base many of the parameter valueson information about ligands and cells in the biological literature. As in many current computa-tional models, these values come from diverse species and experimental conditions, and should beinterpreted with caution.Table 6.1 lists the model parameters, their units and values. (See Supporting Information A.3.3for a discussion on parameter estimation.)6.3.4 Discrete cell modelThe detailed implementation of the HyDiCell3D is described in chapter 3. Briefly, cells experiencemechanical forces from (1) mutual interactions (very-short-range repulsion for excluded volumeeffects and short-range adhesion to neighboring cells; both these forces depend on the distancebetween neighbor surfaces along a line segment joining neighbor cell centers, (2) chemotaxis upgradients of CXCL12a and FGF, (3) drag forces, and (4) adhesion to the substrate (normal to thedirection of motion, and hence only affecting the flatness of the PLLP, but not its motion). Thebalance of these forces determines cell shape and movement. Model cells are both sources andsinks of Wnt and FGF ligands (tracked using reaction-diffusion equations solved on a 3D Cartesiangrid above the plane of the substrate). The slow diffusion of the CXCL12a chemokine in the planeof the substrate is tracked in 2D. I assume that secretion and uptake of ligands in a small volumeelement (voxel) is proportional to local cell-surface area in the given voxel.The receptor levels WR,FR of each cell are computed according to Eqs. (6.4). Those cells withWR >Wthres are then assigned to be a Wnt receptor (WntR) expressing cell (shown red in Fig. 6.3)versus cells with FR > Fthresh denoted FGF receptor (FGFR) expressing cells (colored green). Ifboth of these conditions hold then the cells are both FGFR and WntR expressing cells (coloredyellow) neither of these conditions hold, the cells are taken to be of undetermined fate (colored60Table 6.1: Summary of model parameters.Parameter Definition Units Value ReferenceF0 FGF IC50 nM 0.1 setF1 steady-state FGF receptor concentration nM 1 setW0 Wnt IC50 nM 0.001 setW1 steady-state Wnt receptor concentration nM 0.03 [48]n = m exponents dimensionless 5 setδF FGF ligand decay rate min−1 0.002 [5]δW Wnt ligand decay rate min−1 0.03 [45]rF FGF receptor endocytosis rate min−1 1 setrW Wnt receptor endocytosis rate min−1 1 setkF,on FGF ligand binding rate nM−1min−1 0.1 [65]kW,on Wnt ligand binding rate nM−1min−1 0.005 [48]DF FGF diffusion coefficient µm2 min−1 600 estimatedDW Wnt diffusion coefficient µm2 min−1 600 estimatedpF FGF secretion rate for W1 gradient min−1 24 setpW Wnt secretion rate for W1 gradient min−1 4 setb W1 gradient slope nM µm−1 0.01 setpF FGF secretion rate for WR gradient min−1 15 setpW Wnt secretion rate for WR gradient min−1 5 setb WR gradient slope nM µm−1 0.1 setgrey). The threshold values Wthres,Fthres are assumed to be 90% of the steady state receptor levels(denoted W1,F1, see Table 6.1). The results in Fig. 6.3 are from simulations using the full 3Dmodel but at an early time point (30 minutes of simulation time) when the PLLP has polarized butmigration has not started.Wnt (respectively FGF) expressing cells act as sinks that bind the given ligand. I assume that allWnt receptor expressing cells secrete FGF and Wnt ligands. To compare the results to experimentalobservations, I also show the level of bound receptors, WB and FB on the simulated cells, in Fig. 6.3right panel. WntR active cells, WB, are shown in pink, FGFR active cells, FB, cells are shown inpurple and cells that are neither are shown in grey.A stripe of CXCL12a is produced from head to tail by underlying muscle cells, undergoes slowdiffusion in the plane of the substrate, and is degraded by PLLP cells as they migrate above it.Its concentration is simulated along the 2D surface on which the PLLP moves, according to theequation:∂C∂ t= DC∇2C+α−δCC−κCWCΩW −κCFCΩF (6.6)where DC is the diffusion rate of CXCL12a, α is a constant production of CXCL12a along thecentre line of the domain, δC is the uniform nonspecific degradation of CXCL12a, κCW (κCF ) isthe breakdown of CXCL12a by cells expressing WntR (FGFR/undetermined). ΩW (ΩF ) represents61Figure 6.3: Illustration of the simulations using the 3D deformable ellipsoid cell-basedmodel. View from above of a 1-cell layer representing the PLLP. The stripe of CXCL12awould be directly underneath the PLLP, not explicitly shown. For visualization purposes,the Wnt (blue) and FGF (yellow) ligands are shown as a pair of clouds in each split-image. Left panel: cells are colored based on receptor expression level (red for high Wntand green for high FGF receptor levels, WR,FR). FGFR and WntR expressing cells arecolored yellow. Right panel: cells colored by their bound receptor levels (pink for highWnt and purple for high FGF bound-receptor levels, WB,FB). In the model, the latter is tobe interpreted as the Wnt or FGF signalling levels. Cells in the back of the PLLP expressFGF receptors but do not signal since FGF ligand is so low that most FGF receptorsare unbound. Grey or yellow cells at the back of the PLLP are those that are not yetcommitted to being either WntR or FGFR active cells. Results from the full 3D modelafter 30 min of simulation time (before the onset of migration).the local fraction of cells expressing WntR (FGFR/undetermined), with 0 ≤ ΩW,F ≤ 1. I chose thediffusion, decay and production rates such that in the absence of the primordium there is a steadystate level of CXCL12a along the centre line of the domain.Cells expressing FGFR and WntR respond to gradients of FGF and/or CXCL12a, provided theconcentrations of the ligand(s) are above some detectable level. In that case, the cell orients towardsa vector sum of the gradients, and moves with constant speed (chemotactic force is assumed to havea constant magnitude, independent of ligand type and gradient steepness). If ligand is undetected,the cell’s axis of polarity (and hence direction of motion) can reorient randomly by up to 20◦ in atime step of 1 min (see Table 6.2, for parameter values).Overall simulation protocolIn all cases, I initiate the system with a rectangular configuration of 100 uncommitted (grey) cells,and no ligand of either type. I first simulate the system with cell-cell, cell-CXCL12a adhesion andcell contact forces, but no FGF nor Wnt ligand. This allows the transients due to initialization to beresolved. The time-stepping proceeds as follows:62Table 6.2: Summary of 3D Model Parameters.Parameter Definition Units ValueµECM Drag dyne per s/mm 10−3r Cell radius µm 5Fadhesion Cell-Cell adhesion nN 0.8Factive Force of chemotaxis nN 0.7Fexclusion Exclusion force nN 50Random force nN 0.1Cell-surface adhesion nN 1.1DC CXCL12a diffusion µm2 min−1 20α CXCL12a production nM 0.2δC CXCL12a external decay min−1 0.2kCW CXCL12a decay by WntR cells min−1 0.05kCF CXCL12a decay by FGFR cells min−1 1.0pW Wnt secretion nM 3.0pF FGF secretion nM 0.5FGF external decay min−1 0.15Wnt external decay min−1 0.15KF FGF M-M constant nM 1.0KW Wnt M-M constant nM 1.01. The local concentration of the signalling molecules around each cell is calculated.2. The secreted ligands are distributed into the lattice cubes and the diffusion and depletion ofFGF, Wnt and CXCL12a is calculated (see Eqn. (A.1)).3. Cells orient towards the FGF gradient (green, FGF expressing cells) or CXCL12a gradient(red, Wnt expressing cells) if ligand is at a detectable level. Otherwise, cells orient in arandom direction, biased towards the previous direction of motion.4. All the active and exclusive forces acting on each cell are calculated.5. Cells are moved according to the equation of motion (Eqn. (3.3)).6. This is repeated for every time step.6.4 Tissue subdivision into leading and trailing zonesThe initial step was to explore mechanisms for the subdivision of the tissue into a Wnt-signallingleading domain and a FGF-signalling trailing domain. A specific goal was to determine how aspatially graded parameter, such as the steady state level of the Wnt receptors (W1), across thelength of the primordium can set up this domain subdivision, and where the domain boundary63would be positioned. First, I explore this from a mathematical perspective, linking the outcome tounderlying biochemistry, and then show the results numerically. A 1D model simplification is usedfirst, followed by a full 3D simulation for the PLLP.6.4.1 Theoretical considerationsThe separation of timescales of ligand-receptor binding (fast) and receptor presentation/endocytosis(slow), allows the use of Michaelis-Menten kinetics to simplify Eqs. (6.4) so that the model is recastin terms of the total level of cell-surface receptors, FR(x, t),WR(x, t) and the ligand concentrationsF(x, t),W (x, t).In absence of the negative feedback, the steady state receptor level is (FR,WR)= (EF,1/rF ,EW,1/rw)≡(F1,W1). The variables are then scaled by these steady state levels. Under this scaling, Eqs. (6.4)can be rewritten in the formdFRdt= rF[11+(WR/ω)m−FR](6.7a)dWRdt= rW[11+(FR/φ)n−WR](6.7b)whereφ =F0F1(KF +FF), ω =W0W1(KW +WW)(6.7c)(see additional calculations in Supporting Information A.3.1, and note that FB(x, t), WB(x, t) havebeen eliminated from the equations). In Eqs. (6.7), the space and time dependent quantities φ(x, t),ω(x, t) are related to underlying molecular biochemistry such as ligand concentrations, bindingconstants, and signalling IC50 parameters.The behavior of Eqs. (6.7) can be understood by sketching nullclines of the system (curvesthat depict loci for dFR/dt = 0, and dWR/dt = 0) as shown in Fig. 6.4. Steady states are pointsof intersection of those nullclines. This is done most easily in the “sharp switch” limit of theseequations where Hill functions are approximated as on-off switches (n,m→ ∞), as the analysis ofsteady states is particularly simple [73]. The steady states will be at intersections of piece-wiseconstant nullclines in this limit, namely,FR nullcline: FR = 1 if WR < ω; FR = 0 otherwise,WR nullcline: WR = 1 if FR < φ ; WR = 0 otherwise.The spectrum of behavior of Eqs. (6.7) can be analyzed for constant values of the quantitiesφ and ω . Subsequently, the behaviour of the system when these quantities are time and spacedependent can be described. In particular, as concentrations and other biochemical parameters64Figure 6.4: Analysis of mutual inhibition dynamics leads to model insights. The scaledmutual inhibition Eqn. (6.7) can have one of four possible behaviours represented bynullclines in the WRFR phase plane. (FR nullcline in black, WR nullcline in red.) (a) Bista-bility, in which either a high Wnt receptor (WR) or high FGF (FR) receptor expressingstate can result, depending on initial conditions. (b) A Wnt-receptor expressing state al-ways result, (c) a coexistence-state, with both Wnt and FGF receptors expressed, and (d)an FGF only receptor expression state exists. When the Hill coefficients n,m are large,the steady states (appropriately scaled) occur approximately at a subset of the points{(1,0),(1,1),(0,1)} and transitions between the four qualitative outcomes displayed inthis figure can be summarized by simple inequalities in terms of aggregate quantities φand ω , see Eqn. (6.7c). Figure courtesy of Leah Edelstein-Keshet.vary across the primordium, both φ and ω will also vary across the domain. In principle, severalpossibilities are theoretically possible for the mutual-inhibition model in the spatial context. Forexample, opposed gradients of φ and ω could lead to a transition from Wnt-dominated front, toFGF-dominated rear with or without an overlap region in the middle [82]. This would occur if eitherφ or ω cross their respective thresholds (φ = 1, ω = 1) somewhere in the domain (a coexistencestate, which implies a smooth transition zone between front and back requires that φ > 1,ω > 1over some region in the middle of the domain). Within the context of the simplified model for thePLLP derived in this research, this scenario cannot happen because opposite gradients of φ and ωdo not occur. (This stems form the fact that both Wnt and FGF ligands are graded from front toback.) As will be discussed in the next section, the domain always evolves to the bistable state(Fig. 6.4a) with φ < 1,ω < 1 everywhere. The front becomes WR dominated, and the rear becomes65FR dominated, but the transition zone is sharp, and determined by a separatrix (boundary of twobasins of attraction). Future models that consider the role of Dkk1 and/or other effectors that varyspatially might amend this feature of the model.6.4.2 The 1D spatial simulationThe signalling model (Eqs. (6.5) to (6.7)) is implemented in one spatial dimension, with defaultparameters estimated from the literature, where possible (Table 6.1). First, the length of the PLLPis scaled to L = 1, with rear and front at x = 0 and x = 1, respectively. To understand how thesignalling domains might form, initial conditions in the coexistence state are considered (Fig. 6.4)with W (x,0) = 0.01, F(x,0) = 0.005, WR(x,0) = 0.1, and FR(x,0) = 0.01. These represent smallinitial values of the scaled ligand and receptor levels (scaling previously discussed) that would tendtowards the unique coexistence state over time. The parameter W1 is assumed to initially be gradedacross the PLLP (W1(x)= 0.01x+0.03). Motivating this assumption is the experimental observationthat the Wnt signalling dependent transcription factor lef1, has graded expression (highest at front).Next, the rate of secretion of Wnt and FGF ligands are adjusted to obtain biologically realisticsignalling domains. Signalling domains form as expected. Note that zones readily form with manyother types of initial bias, e.g., graded WR, see Supporting Information in [47], or flux of Wnt ligandfrom the right boundary.Results of the 1D simulations are displayed in Fig. 6.5. Kymographs of the ligand and receptorlevels of each type (Wnt, top, and FGF, bottom) across the 1D length of the tissue, are shown.After some transient behaviour that depends on initial conditions, the formation of robust zones isobserved by roughly 20 minutes, with strong interfaces that persist. The shapes of the chemicalprofiles and receptor densities at t = 25 minutes are shown on the left panel of Fig. 6.6. BothWnt and FGF ligands are graded towards the front, with FGF ligand the sharper of the two. Thereceptor densities exhibit a sharp interface in the middle of the domain, as expected from the mutual-inhibition model.An investigation of the time evolution of the values of φ and ω over the domain was performed.At t = 0, given the initial gradient in W1, every point in the domain satisfies φ > 1,ω > 1, implyinga single coexistence (Wnt and FGF) steady state everywhere. The gradient in W1 induces a gradientin ω , which ensures that less Wnt activity at the front of the PLLP is required to inhibit FGF activity.As FGF and Wnt ligand is produced, ω and φ fall below 1, while WR and FR signalling levels evolve.First, as ω decreases below 1, the PLLP transforms to a Wnt-dominated system. Second, as φ dropsbelow 1, the PLLP becomes a bistable system, with spatially varying FR and WR conditions. Whilethe front is already locked into Wnt dominance, the rear is diverted from its Wnt-wards trajectoryinto a basin of attraction of FGF-dominance, and becomes FR dominated. In the final state, a sharpzone boundary forms between front and rear, determined by parameters that define the separatrix inthe bistable state.66Figure 6.5: The Wnt-FGF mutual inhibition model predicts the formation of signallingdomains in response to graded W1 across the primordium. Simulation kymographsfor Wnt and FGF signalling. Initially, a gradient W1(x) = bx+0.03 is assumed in the un-coupled Wnt steady state parameter, with b = 0.01. Time increases downwards; positionacross the PLLP is horizontal, with leading edge on the right. (a) Wnt ligand, W (x, t),(b) Wnt receptors, WR(x, t), (c) bound Wnt receptors, WB(x, t), (d) FGF ligand, F(x, t),(e) FGF receptors, FR(x, t), and (f) bound FGF receptors, FB(x, t). Signalling domainsform after a few minutes, with a sharp boundary between zones. Bound ligand concen-trations are calculated using Eqs. (6.3). Parameters are as in Table 6.1. Initial conditions:W (x,0) = 0.01, F(x,0) = 0.005, WR(x,0) = 0.1, FR(x,0) = 0.01. Figure courtesy of ColeZmurchok.6.4.3 Parameter dependence of the zone boundaryNext, the effect of the various parameters in the signalling model on the position of the zone bound-ary was investigated. Since the competition between Wnt and FGF signalling is biased by the ligandconcentration, by ligand binding to receptors, and by the sensitivity of a given cell type to boundreceptors of the opposite type, it stands to reason that any one of these parameters would affect thesubdivision of the tissue into zones.To investigate this dependence, the 1D simulations were conducted, as before, and the locationof the sharp receptor interfaces for 20 values of each of the parameters W0, F0 (“IC50” cross inhibi-tion sensitivities), F1 (steady state FGF receptor level), pF , pW (ligand secretion rates), and b (slopeof initial W1 gradient) were identified. The boundary location is shown in Fig. 6.6 as a functionof each of these parameters. The top edge of the panel represents the leading edge (x = 1), while67Figure 6.6: The position of the boundary between leading and trailing zones is parameterdependent. In the left panel, the concentration profiles of the ligand and receptor activityare plotted at t = 25 minutes (after formation of signalling domains), which reveals theformation of a sharp boundary between the leading and trailing zones. In the right panel,the boundary position between the leading and trailing zones is shown to be a functionof model parameters. Parameter sweeps were conducted for a range of “IC50” crossinhibition sensitivity parameters F0,W0, FGF steady state receptor levels F1, productionof Wnt and FGF ligand pW , pF , as well as the slope of the W1 gradient across the PLLP,b. The vertical axis represents the fraction of the PLLP length that is FGF-receptorexpressing at steady state. The horizontal axis represents the given parameters, scaled bytheir baseline values as listed in Table 6.1. Figure courtesy of Cole Zmurchok.the bottom (x = 0) is the trailing edge. When the boundary position is 0, the entire primordium isdominated by Wnt receptor expression, and when the boundary position is 1, the entire primordiumis dominated by FGF receptor expression. The horizontal axis gives the ratio of parameter value tothe baseline value used to produce signalling domains as in the left panel of Fig. 6.6. For example,to produce the signalling domains in the left panel of Fig. 6.6, b= 0.01, and the boundary position isapproximately at x = 0.5. When this parameter is increased to b = 0.02, the ratio between b and itsbaseline value is 2, and the boundary position can be identified from the light blue curve in Fig. 6.6to be approximately at x = 0.25.As expected, increasing either the basal (uninhibited) steady-state FGF-receptor level, F1, or the“IC50” inhibition parameter W0 (which decreases the sensitivity of FGF receptor presentation to in-hibition by Wnt) biases the cells towards FGF receptors and hence pushes the boundary towards thefront while increasing F0 (which decreases the sensitivity of Wnt receptor presentation to inhibitionby FGF) has the opposite effect. The effects of changing each of the aggregate parameters, W0, F0or F1, on the boundary position between the leading and trailing zones matches my intuition aboutthe signalling dynamics. For example, by increasing F1 or W0, the PLLP is biased towards FGFreceptor expression. If F1 increases, the parameter φ decreases in Eqs. (6.7). This reduces the total68number of FGF receptors required to inhibit Wnt receptor expression, whence the PLLP is biasedtoward FGF receptor expression. Similarly, if W0 increases, so does the parameter ω . This increasesthe total number of Wnt receptors required to inhibit FGF receptor expression, whence the PLLP isbiased toward FGF receptor expression. The secretion of Wnt and FGF ligand also influences theposition of the boundary. W and F appear in ω and φ . Hence any changes in the secretion rate pWand pF will affect the number of signalling receptors required to inhibit receptors of the oppositetype. Lastly, note that the steepness of the gradient of W1 (basal steady state level of Wnt receptorswhen there is no FGF inhibition) across the PLLP can affect the boundary position. One observationis that as b→ 0, W1(x)→ 0.03. In this case, no species can induce a gradient in φ or ω across thePLLP, whence no signalling domains form.6.5 Simulation results using HyDiCell3D6.5.1 Formation of zones in 3D simulationsHaving obtained zone partitioning in the 1D model, I next asked whether a similar subdivisioninto leading and trailing domains would take place in the full HyDiCell3D. Consequently, I ran the3D model with the basic parameter set. In all cases, the system was initiated with 100 uncommittedcells, and no ligand of either type. To resolve transients in the initial shape, the cells are first allowedto adhere and interact with each other and with the center of the domain where the CXCL12a stripeis located. This leads to an elongated cluster. (In the absence of the added adhesion to the CXCL12astripe, the cell cluster is round). I then turned on the signalling dynamics with W1 graded acrossthe length of the primordium (W1 = 0.03+ 0.006x, where here, x represent the position within thePLLP in units of µm.) The characteristics of the cells (whether WntR expressing, FGFR expressingor uncommitted) are determined by the levels of FGF and Wnt receptors; an undetermined stateoccurs when both Wnt and FGF receptors are below a set threshold.Results are shown on the left side of Fig. 6.3 (with green for FGFR expressing cells, red forWntR expressing cells, yellow for both WntR and FGFR expressing cells and grey for uncommittedcells). (In all figures, “front” is on the right, “rear” on the left). I also show the Wnt ligand con-centration profile (blue) and FGF ligand (yellow); both ligands are actually distributed throughoutthe 2D myoseptum on which the PLLP migrates (the z = 0 plane), but for visualization purposes,each one is shown as a split image. The right side of Fig. 6.3 is from the same simulation as the leftside but I used different colors to emphasize the distinction between receptor expressing cells andsignalling cells/bound receptors. The pink cells on the right are the WB cells, i.e. cells with highlevels of Wnt ligand-bound Wnt receptor. The purple cells on the right are FB cells, i.e. cells withhigh levels of bound FGF receptors. The grey cells on the right express FGF receptors but sincethe FGF ligand level, and hence bound FGF receptors, are very low, those cells are not consideredFGF signalling cells. These results can be compared to experimental results where lef1 is used as a69marker for the Wnt signalling zone, pea3 is used as a marker for the FGF signalling zone. Wnt10ais one of the candidates for the Wnt ligand contributing to this signalling network, FGF10a is oneof the FGF ligands and the FGF receptors can also be labelled, see Fig. 6.7.wnt10alef1afgf10afgfR1pea3Figure 6.7: Experimentally observed Wnt and FGF expression and signalling levels. RNAin-situ hybridization showing expression at 32 hours post fertilization in the PLLP. Fromtop to bottom these are: wnt10a (Wnt ligand), lef1 (Wnt signalling), fgf10a (FGF ligand),FGFR1 (FGF receptors) and pea3 (FGF signalling). The scale bar is 10µm. Figurecourtesy of Damian Dalle Nogare.I ran several simulations in which the aggregate parameters F0,W0,F1,W1 were varied by a factorof 10 up and down from their baseline value (represented with an asterisk, X∗, see values in Table6.1) to study the effect on the resulting domain subdivision. Results after 50 min of simulation time(prior to migration) are shown in Fig. 6.8. Trends match Fig. 6.6, indicating that the analysis of the1D model parallels the observations for the 3D model. The ligand profiles for both FGF and Wntare also in good agreement with results from the 1D model.70Figure 6.8: The position of the boundary between the leading and trailing zones is pa-rameter dependent. The four parameters, F0, F1, W0 and W1 are altered. In the firstcolumn the baseline value (X∗) has been multiplied by 0.1 and in the last column it hasbeen multiplied by 10. Increasing F0 (top row) and W1 (bottom row) results in an increasein the size of the Wnt domain. Conversely, increasing F1 (second row) and W0 (third row)results in a decrease in the size of the Wnt domain. Results from the full 3D model after50 min of simulation time (before the onset of migration).6.5.2 MigrationHaving identified conditions for the formation of leading and trailing zones, I next asked whatfurther conditions are required for collective migration of the PLLP. To probe this question, I as-sembled basic facts from the biology (surveyed above) supplemented, where needed, by reasonableminimal assumptions. In actual PLLP experiments, it is known that aside from the Wnt-FGF polar-ization of the tissue, the front and rear portions are also distinct in levels of chemokine receptors.Cells in the leading 60% of the PLLP express high levels of CXCR4b, whereas those in the remain-ing rear portion express high CXCR7b, [2]. This does not exactly coincide with the growth-factorzones: the leading 1/3 of the PLLP is the WntR activity zone [3]. However, here the simplification71is made that from the standpoint of interactions with CXCL12a, the WntR-FGFR subdivision servesas reasonable proxy for the chemokine zone subdivision. One reason for this simplifying assump-tion is that I thereby avoid the further complications of untangling the signalling circuits that set upthe zones of differential chemokine receptor expression. A second reason is that little is to be gainedby this intermediate step currently, when the detailed biology of this step is largely unknown.Basic migration modelTo model the migration of the cells (Fig. 6.9), I assumed a narrow long stripe of CXCL12a thatwas continually produced along the central long axis of the domain, as observed in experiments[17, 38, 95]. I further assumed the following: (a) Both leading and trailing cell types degradeCXCL12a, but trailing (FGFR expressing, green) cells do so at a 20-fold higher rate than leading(WntR expressing, red) cells ([12], and Fig. 1 in [52]). (b) Leading cells are assumed to migrate bychemotaxis up CXCL12a gradients. (c) FGFR expressing cells (shown in green) are pulled by theiradhering neighbors and move by chemotaxis up FGF gradients. I consider alternative hypotheses inthe Discussion, section 6.7.Results of the basic simulation (Fig. 6.9) are shown at four time points, in the very beginningof the simulation, after 5, 100 and 200 minutes. The top panels show the evolution of the PLLPshape, the secretion of both Wnt and FGF ligands, the polarization, and the migration across thedomain corresponding to about 250 µm. The speed of motion is the same order of magnitude as theexperimental observations (1-2 µm/min). In the same figure, I also show the profiles of Wnt ligand(blue dashed line), FGF ligand (green dotted line), and CXCL12a (red solid line) across the lateralaxis of the PLLP. It is evident that the Wnt and FGF ligand profiles form a sharp peak close to thefront. I also observe a trench in the CXCL12a profile due to uptake and degradation of CXCL12aby the trailing cells. The simulated CXCL12a profile agrees well with results of [97].Absence of FGF chemotaxis leads to migration failureIn my model, front cells migrate up CXCL12a gradients, and exert pulling forces on their neigh-bours, and cells further back actively migrate up a gradient of FGF ligand. I asked whether activemigration of the rear is essential for collective migration, or whether the cells at the front might be“strong enough” to account for the PLLP motion.I first estimated the viscous drag of a cell crawling through tissue. To do so, I used the factthat in this regime, cell speed is set by a balance between the active pulling forces exerted by a celland the drag forces it experiences (F ≈ µv, where µ is a drag coefficient). In cell collectives andepithelia, typical force generated by a cell has been measured to be on the order of a few nN [14, 22].(This is reasonable, and represents roughly the combined effect of around 1000 actin filaments, eachexerting a pushing force of about 1 pN at the cell edge [29].) The velocity of the cells in the PLLPis observed to be around 1 µm per min [38]. This leads to an estimate of µ ≈ 6 · 10−3 dyne s/µm72Figure 6.9: Successful migration in the discrete cell model. Simulations of the migratingprimordium showing WntR expressing cells (red), FGFR expressing cells (green), WntRand FGFR expressing cells (yellow) and indeterminate cells (grey). Top left: at 2 min,top right: at 5 min, bottom left: at 100 min and bottom right at 200 min. The line graphsshow concentrations of CXCL12a (red), FGF (green) and Wnt (blue). Parameters as inTable 6.1.for the drag coefficient (1 dyne= 10−5 N) and is in the same range as the estimate in [70] for a cellin a Dictyostelium discoideum slug. Cell motion is dominated by viscous drag, so inertial forces areinsignificant and acceleration terms are absent in my model (contrasting with the model in [20]).Using the above reasonable estimate of µ , I found that if only the cells at the front provideforward propulsion, then the PLLP hardly moves at all. Without an additional active propulsionmechanism for trailing cells, such as the FGF-directed chemotaxis, migration was either absent, ormarginal at best. Furthermore, the PLLP became unrealistically elongated and the front occasionallytore apart from the rear (see Fig. 6.10). At the same time, if cells in the rear exert overly large73motile forces relative to those at the front, then the PLLP shape becomes too wide and distorted,in disagreement with experimental observations. This points to the need for proper balance of themotile forces in the front and back of the PLLP.Figure 6.10: Abolished FGF chemotaxis. If FGF expressing cells in the trailing part of thePLLP do not chemotax towards FGF ligand, then the PLLP becomes elongated, and itsshape no longer resembles the experimental observed PLLP shape.6.5.3 Cell growth and proliferationIt is well known that cell growth and division take place during the course of PLLP migration withthe total number of cells doubling throughout the process [3, 94]. This fact led me to inquire howsuch cell proliferation contributes to the zone formation, migration, and overall shape of the PLLP.While it was originally thought that leading zone cells proliferate faster [3, 94], more recent exper-iments revealed a relatively constant and spatially uniform cell cycle length (539 ± 127 minutes)across the PLLP [57]. In the simulations, the cell’s rate of growth in each time step is chosen from auniform distribution between 0 and 0.001 min−1 (corresponding to a cell cycle between 0 and 1,000minutes). Cell division occurs when the volume of a given cell has doubled. The daughter cellshave the same parameter values as the mother cell.Results with a range of growth rate values are shown in Fig. 6.11. For the higher growth rate,I found that as the PLLP migrates and elongates, the cells at the back no longer detect FGF, andeventually break away from the rest of the PLLP.6.6 Comparison with published experimentsThe spatiotemporal FGF and Wnt ligand distribution is hard to quantify experimentally. Hence, di-rect comparisons of model predictions to experimental observations are not possible for the chemicaldistributions of FGF and Wnt ligands. For this reason, I sought simple experimental comparisonsthat could (in)validate the computational model. I consequently compared the model behaviourto the results of laser ablation (cutting) experiments, even though such experiments have alreadybeen successfully recapitulated in a computational model [16]. I further simulated several mutants74Figure 6.11: Effect of growth rate. Simulation results for three different growth parameters.Simulation time is 6 hours. As the growth rate is increased, the shape and size of thePLLP becomes distorted, and fails to match experimental observations.to check whether the model predictions agree with what has been observed. In the experiments,staining has been done for Wnt signalling cells, using lef1 as a marker, and for FGF signalling cells,using pea3 as a marker. Hence, the cells in the simulations are coloured based on signalling cellbehaviour.6.6.1 Laser ablation in silicoTo compare the model to results of [16], I performed ablation (cutting) experiments in which part ofthe simulated PLLP was removed. Results are shown in Fig. 6.12. The left-most column shows thecontrol PLLP before ablation, the second column shows the configuration at 12 min, directly afterablation, and the next two columns show what happens to the PLLP as a result.I performed several such in silico experiments. Following the establishment of clearly definedleading (pink) and trailing (purple) zones of Wnt and FGF activity, respectively, removal of theleading Wnt active zone (Fig. 6.12, Row 1), initially leaves behind a PLLP entirely composed ofFGFR active cells. However, in the absence of fresh FGF or Wnt signals produced by leading cellsthe entire remaining PLLP assumes an undefined state (grey) characterized by lack of both WntRand FGFR activity. Unsurprisingly, the same behaviour occurs when I remove both front and back,see Fig. 6.13 . Removing the trailing cells, but leaving the middle fragment that contains FGFRactive cells (Fig. 6.12, Row 2), enables the PLLP to continue migrating. This result is comparable75Figure 6.12: Recovery after partial ablation. As in laser ablation from [16], part of thePLLP is removed in the simulation, and recovery observed. Top row: the front segment(WntR active cells) is removed, leaving only FGFR active cells; the PLLP stalls andcannot continue. Middle row: the rear portion is removed, leaving a few FGFR cellsbehind; the PLLP migration continues. Bottom row: the middle segment is removed.Motion is stalled while the trailing FGFR active cluster catches up with the front. Oncethe clusters have merged, migration resumes.to the experiments where staining shows that the middle fragment has FGFR activity and the PLLPmigration stalls temporarily and then starts up again. If I remove a larger part of the trailing PLLP,leaving no FGFR active cells behind the PLLP has no net forward movement, see Fig. 6.13. Aboutthe parameter W1 I assumed the following: Initially, W1 is graded across the primordium, but whencells start dividing, daughter cells inherit the value of W1 from their mother cell.When the middle portion of the PLLP is removed (Fig. 6.12, Row 3), the back catches up tothe front by chemotaxis towards the FGF gradient created by leading cells. The front and rearclusters merge, creating a gradient of CXCL12a and migration continues. These predictions are inqualitative agreement with both experiment and computational model in [16].6.6.2 Comparison with Mutants and other experimentsI compared the results of the simulations with experiments described in [2]. They used heat shock tooverexpress the Dkk1 inhibitor and found decreased expression of FGF3 and FGF10 ligand. Stain-ing shows no FGF activity (using pea3 as a marker) and no Wnt activity (using sef as a marker).76Figure 6.13: Laser ablation experiment Top panel: Both the front and the back of the PLLPare removed resulting in the middle segment stalling. Bottom panel: A large part of theback is removed, leaving no FGFR active cells behind. The PLLP has no net forwardmovement.Inhibiting Dkk1 experimentally is comparable to decreasing F0 in the simulations, the IC50 param-eter for Wnt inhibition. When F0 is lowered, fewer bound FGF receptors (and hence also less FGFligand) is needed to inhibit Wnt receptor presentation Eqn. (6.4c). The model results (Fig. 6.8)show that decreasing F0 10 fold from the baseline level suppresses Wnt and FGF activity: there areno WntR active (pink) cells and no FGF ligand (yellow cloud) as a result.The receptor blocker drug SU5402 reduces the ability of the FGF receptor to respond to FGFligand. According to experiments in [2], after treatment with the SU5402 inhibitor, WntR signallingactivity takes over the entire primordium. Furthermore, using the SU5402 inhibitor just prior tomigration leads to stalling of the primordium. In the model, this treatment corresponds to decreasingthe steady state level of FGF receptors, F1. As shown in the third row of Fig. 6.8, decreasing F1increases the WntR activity domain and stalls migration, in agreement with experiments in [2].Wnt ligand interacts with its receptor and inhibits APC-dependent β-catenin degradation to pro-mote Wnt/β-catenin signalling. The apcmcr mutation reduces degradation, stabilizes β-catenin andpromotes Wnt/β-catenin signalling independent of the Wnt ligand. As a result, FGF signalling isunable to inhibit Wnt signalling by promoting expression of a diffusible inhibitor of Wnt, Dkk1. Un-regulated Wnt signalling in apcmcr mutants promotes high levels of Wnt-dependent FGF expressionin the PLLP. Wnt signalling is unable to inhibit the FGF signalling that results from this unusuallyhigh levels of FGF expression, and the PLLP is left in a state characterized by high levels of bothWnt and FGF signalling, where mechanisms that normally determine FGF-Wnt mutual inhibitionare lost. To reproduce these results in the simulations, I turned off FGF-Wnt mutual inhibition,which decouples the biochemistry of neighbouring cells. Thus, as the 1D model predicts, all the77cells express high levels of FGFR and WntR activities.In addition to these previously published results, new dose-response experiments were con-ducted by my collaborator (Damian Dalle Nogare, NIH) to compare to simulation results of Figs. 6.6and 6.8, as described in Supporting Information A.3.2. Using the FGF receptor blocker SU5402 heshowed that the size of the WntR activity domain increases with increasing amount of the inhibitor,Fig. 6.14 (in good agreement with similar experiments for a different cell type [93]). These resultsvalidate the predictions of the model where decreasing the parameter F1 (the steady state level ofFGF receptors in the absence of inhibition) increases the size of the WntR activity zone.Figure 6.14: Dose response experiment using SU5402. The FGF receptor blocker SU5402increases the size of the WntR activity domain in a dose-dependent manner (lef1 usedas a marker). The size of the lef1 domain has been normalized relative to the size of thewhole PLLP. A Linear regression analysis resulted in y = 0.2005x+0.515. A one-wayANOVA reveals that the slope is significantly non-zero (0.2005 ± 0.03012). The 95%confidence interval for the slope is 0.1405 - 0.2605, with R2 = 0.3536, and P< 0.0001.These results validate the findings of my model. Figure courtesy of Damian DalleNogare, NIH.6.7 DiscussionHere, I investigated the lateral line primordium development and migration in zebrafish. In thisintriguing system, a cluster of 100-200 cells maintains cohesion and directionality while migratingover a distance one order of magnitude longer than the cluster length. First I showed that the78mutual inhibition of Wnt and FGF signalling was able to, on its own, set up the initial polarizationof the PLLP into distinct leading and trailing domains under appropriate assumptions (researchquestion 1). This model links ligand-receptor binding (in the Michaelis-Menten regime) to theinhibitory effect (of bound receptors on the antagonist receptors). Thereby, predicting how theboundary between Wnt and FGF domains depends on parameters such as steady state receptor levels(in absence of inhibition) and sensitivity to inhibition by the antagonist bound receptors (researchquestion 2).The theoretical model discussed here has more universal applicability, and variations on thistheme could be at play in other differentiation processes of tissue subdivision. In particular, “crossedgradients” (i.e. gradients of opposite slopes) in the two ratios of receptor-ligand properties, iden-tified as φ and ω , could produce a tissue with a central zone of overlapped influence (transitionacross the bottom panels of Fig. 6.4), or with a sharp boundary (transition across the top panels ofthis figure). For this model of the PLLP specifically, such crossed gradients do not exist because thespatial variation of φ and ω stem from the FGF and Wnt ligand distributions, which are both gradedin the same direction (increasing towards the front). Nevertheless, tissue subdivision occurs, oncethe PLLP evolves into a bistable state, with WntR activity locked into the front, and FGFR activityovertaking WntR in the back.I showed that the CXCL12a adhesive stripe provides directionality, as well as, coupled to activedegradation by trailing cells, forms the guidance cue that directs PLLP migration along the lateralaxis of the embryo. Given the viscous drag, it appears highly probable that all cells in the PLLPare actively motile. Simply having a “front engine” (a few cells at the front) pulling the entire masswould lead to migration stalling, and even fragmentation of the group, which is not observed undernormal conditions (research question 3). Chemotaxis is the simplest and most plausible mechanismfor the active motion of trailing cells, though not excluding the possibility that trailing cells alsoactively follow leader cells through some other mechanism (such as micro-protrusions).The shape and speed of the collective migrating mass depends on numerous factors such as rateof cell growth (which determines cell division rate), active crawling and adhesion of cells to one an-other and to CXCL12a. When these factors are increased or decreased beyond certain ranges, shapeand speed are abnormal, recapitulating a variety of mutants. Wnt signalling regulates cell growthrate, and FGF signalling is implicated in cohesion. Both signalling types determine differentiationand the distinct roles of leading and trailing cells. WntR active front cells sense the CXCL12a gra-dient and steer the cluster. Rear cells eventually express CXCR7b chemokine receptors and degradeCXCL12a. The model demonstrates that this could form the self-propagating gradient to whichWntR cells orient. FGFR active cells eventually engender the initiation of neuromast depositionthat has not yet been addressed in this research.It is interesting to compare the results in this chapter with those of previous modelling papers[1, 12, 16, 20, 88]. In Streichan et al. [88], the PLLP is modelled as a 1D rod, moving through a self-79generated gradient of a single chemokine (CXCL12a), similar in flavor to “treadmilling” solutionsin [23]. The rod binds the diffusible chemokine and moves with velocity assumed to be proportionalto the gradient so created. The model for chemokine profile and rod position sustains a travellingwave solution, representing the steady-state long-range PLLP migration. The authors further con-sidered deposition, assuming that below a threshold gradient of CXCL12a, the rear cells in a discrete1D elastic array fall off. They derive an analytic expression for the migration velocity as a functionof rates of ligand degradation and diffusion, the size of the PLLP and the rate of chemokine uptake.This model provides elegant analytic insights into minimal requirements for “self-generated” gra-dient chemotaxis in a general setting, while simplifying the biology and neglecting key features ofthe PLLP. For example, the authors do not consider interactions between CXCR7b and CXCL12aligand, which is now well established in the literature. In comparison, my model begins to addressthis limitations by including some of these biochemical interaction, consequently limiting many ofthe results to numerical experiments.My collaborators, Ajay Chitnis and Damian Dalle Nogare, at the National Institute of Health,previously crafted a Netlogo model [12, 16] as an aide to decipher laser ablation experimentalobservations. The simulation uses “turtles” attached to one another with springs: when one moves,its neighbor moves with a slight lag, disallowing overlaps. This simulation had the advantage ofclose integration with experiments. While remaining fairly elementary, it was able to capture keyaspects of the biology, and explain experiments. A limitation was the absence of realistic forcesof cell-cell interaction and friction. Consequently, a conclusion, based partly on the assumption ofpulling exerted by few leading cells, was that cells respond to the CXCL12a concentration ratherthan CXCL12a gradient. In the simulation here, I endow a larger number of cells with the abilityto sense, orient to, and migrate up the CXCL12a gradient to correct that feature. Furthermore, asargued previously, drag forces would stall the PLLP were it not for additional chemotaxis motionof the trailing cells.Allena and Maini [1] described chemical signalling coupled to cellular mechanics in the PLLPusing a reaction-diffusion finite element model. Cells are depicted by 2D Maxwell viscoelastic ele-ments inside an elliptical PLLP. The incorporation of additional signalling (FGF and Wnt ligands,CXCR4b and CXCR7b chemokines as “readouts” of the ligand distribution), coupling mechanicsand biochemistry, and the successful recapitulation of mutant behaviours are all strengths of thismodel. As in Allena and Maini [1], I have adopted a hybrid model to describe the PLLP dynamicsin 3D. The discrete cell model builds on limitations of [1] by including interactions between thePLLP and CXCL12a, relative motion and rearrangement of cells, and more direct representation ofintercellular forces. In contrast to [1], my model also explains the initial subdivision of the PLLPinto leading and trailing zones.Another hybrid model [20] is based on cell-cell interaction, cell-substrate adhesion and cell dif-ferentiation. The model consists of discrete cells in 2D with a non-local sensing radius. The PLLP80is described as a “flock” of cells accelerating via haptotaxis to CXCL12a, forces of Cucker-Smalealignment, attraction-repulsion, and chemotaxis to FGF gradient. The model captures the migrationof the PLLP resulting from a graded distribution of CXCL12a. One issue with the representation ofthe PLLP as a “flock” is that macroscopic flocks and swarms generally operate in a high Reynold’snumber regime, whereas I have chosen to neglect inertial forces and acceleration terms and to con-sider a low Reynold’s number regime for the cellular milieu. One interesting feature of [20] is that,by including cell-cell lateral inhibition in the trailing domain (and foci of FGF), their model can de-pict the formation of stable rosette structures within the PLLP (see their Figs. 9 and 10). In contrast,I have not considered the rosette formation in my model but I base the cell motion and interactionson a closer concordance with recent biological results. In particular, I include downstream effectsof Wnt-FGF signalling, i.e., formation of receptor activity domains, chemotaxis of trailing cells toFGF, and degradation of CXCL12a (SDF1a) in the trailing domain.All models described above, including mine, would benefit from further experiments that couldhelp to probe hypotheses or compare with predictions, in particular, the following experimentswould be of primary interest:I Implant additional PLLP cells into the primordium before migration begins. How would thisaffect the size of the Wnt domain and the migration? In this model, all WntR active cells aresources of Wnt ligand. This assumption will contribute to the size of the signalling domains.In the biological system, manipulating the size of the PLLP could validate the role of theWnt10a ligand, or determine if the Wnt/β-catenin pathway remains active at the front of thePLLP.I Similarly, I propose the inclusion of an external source of Wnt ligand (such as a Wnt-ligandcontaining bead) in an experiment with the wild-type primordium. Would this affect the sizeof the Wnt and FGF signalling zones? Can an external Wnt source rescue the signallingdomain and migration in the migration-stalled lef1-mutants?I Repeat the laser ablation experiment in a PLLP mutant with no cellular proliferation. Doesthe removal of the back segment still lead to PLLP migration or does it stall? In the previousexperiments and my model, migration is rescued by the reformation of signalling domains asproliferation sufficiently lengthens the PLLP for the signalling domains to form. This couldaid in the understanding of the role that proliferation and subsequent rosette organization playin maintaining PLLP migration.Primarily, the suggested experiments seek to determine the roles of the Wnt and FGF ligands insustained organization and migration of the PLLP. Does an initial signalling event start the PLLPmigratory process, or, as is assumed in this model, is a continuous source of Wnt ligand required tomaintain the organization and migration of the PLLP?81In summary, this is not the first model of the PLLP migration, though it is the first fully 3Dsimulation of its kind, with realistic cell-cell and cell-substrate forces. While the model is a drasticsimplification of the biology, it nevertheless provides insights about the Wnt-FGF mutual inhibition,forces of propulsion and drag, and accounts for mutants. There are several opportunities for futurework. First, there is the process of neuromast deposition that should be considered in a future study.Second, there are many additional players in the signalling network responsible for segmenting thePLLP into leading and trailing zones, such as Dkk1 and sef, and including these players could aidin the understanding of their specific roles within the PLLP.82Chapter 7The role of adhesion in collectivemigration of cell-sheets7.1 IntroductionThe previous chapter highlighted the complex and important role of collective migration of cellsduring development. In this chapter, I will continue to study collective migration but with a focuson migration of 2D cell-sheets (epithelium) that can be studied in a laboratory setting, in vitro.These experiments provide insights into how cells communicate with each other and with the micro-environment to migrate collectively, an important process in cancer cell invasion and intravasation.I initiated a collaboration with the experimental group of Dr. Calvin Roskelley at UBC toproduce and develop models that are relevant for biologists. This collaborative project is still inits early stages and more experiments and simulations are needed to further our understanding ofepithelium migration. However, I will demonstrate the model-to-experiment feedback cycle thathas been established as well as some preliminary predictions.I will explore the adhesion properties necessary for the emergence of cells at the leading edgeof a cell-sheet that drive the collective migration, called leader cells. The origin of leader cells, aswell as the difference between leader cells and bulk cells, is still largely unclear. Leader cells areobserved at the migrating front of a group of cells and they may appear to drag the group behindthem. They are found to have a front-back polarized biochemical signature and characteristics ofcells that can migrate singly (i.e. mesenchymal-like cells) which includes being spread out withmany long protrusions [58]. It is not clear how leader cells emerge but they have been isolatedexperimentally and shown to have higher activity of certain molecules known to be involved in cellmigration (Rac, β1-integrin and PI3K) than the bulk cells [102].It is not clear how changes in cell-cell adhesion contribute to tumor cell migration through adense cell-sheet towards a leading edge. I will investigate whether a leader-cell phenotype can83emerge in the bulk of the cell-sheet and then migrate to the leading edge. Alternatively, leadercells may be a phenotype that can only emerge in the micro-environmental conditions at the leadingedge. Furthermore, I will ask if changes in cell-cell adhesion molecule expression facilitates theemergence of the leader-cell phenotype when the cells reach the leading edge.I will carry out simulations with a mixture of two cell types with different adhesion properties.The modelling work will guide the design of experiments to be conducted in the Roskelley labora-tory. The experiments are designed to elucidate if differences in adhesion expression can accountfor the leader-cell phenotype taking over the leading edge or being found at the tip of a group ofcells protruding out of the leading edge of a cell-sheet (also referred to as fingers). The classic invitro wound healing assay experiments will be used, where a dense cell-sheet (or monolayer) isgrown on a plate. Subsequently, a scratch or a wound, is made in the sheet with a pipette. The cellsthen migrate into the wound area to close the gap. Wound healing experiments are widely usedand studies with various cell types have shown different morphologies and velocities of the leadingedge.The experiments in the Roskelley lab will be conducted with three cell types that express differ-ent trans-membrane molecules involved in adhesion (cell adhesion molecule (CAM)), called cad-herins. To study the role of different cell-cell adhesion strengths, the experiments will be conductedwith a mixture of two of the cell types. The initial hypothesis is that adhesion bonds between thesame types of cadherins (homotypic adhesion) is stronger than the adhesion bonds between differentcadherins (heterotypic adhesion).Section 7.2 lists the research questions that guide this study. Section 7.3 briefly outlines howthe HyDiCell3D framework is used to set up simulations of cell-sheets that have two cell typeswith different adhesion properties. Section 7.4 shows the results from the simulations. Section 7.5demonstrates the experimental results. In section 7.6, I compare the results from the experimentsto the simulations and show some preliminary results from an algorithm designed to automaticallytrack cells in the experiments.7.2 Research questionsUnderstanding collective migration in 2D sheets and the emergence of leader cells raises somefundamental questions including: Where do leader cells originate? What distinguishes leader cellsfrom trailing cells? How do intra- and inter- cellular signalling affect differentiation and distinctroles of leading and trailing cells? How do cells in a sheet communicate? In my research, I focuson the following key questions1. How do different cell adhesion properties affect the leading edge morphology of an epithe-lium?2. How does proliferation contribute to that morphology?843. How is the rate at which a wound closes affected by cell adhesion properties?4. How do changes in cell adhesion facilitate cell migration toward the leading edge?5. How do the strengths of homotypic and heterotypic cadherin adhesion molecules compare?7.3 Wound healing simulation setupTo study the mechanical attributes of migrating cells in a wound healing assay, I use the HyDiCell3Dframework. In densely packed environments, like the epithelium, the deformable properties ofthe cells allows them to have different numbers of neighbours and the adhesion strength betweenneighbouring cells can vary because the surface area that they share changes. These features havean effect on the migration properties of cells within the sheet.I conduct simulations both with a single cell-type and a mixture of two cell-types with distinctadhesive properties, illustrated with different colors in the simulations. The interfacial force (asdescribed in chapter 3) consists of adhesion between neighbouring cells and an exclusion force torestrict the overlap of cells. Initially the cells are spherical and depending on the forces acting onthem, their shape can become elliptical. The cells all adhere in the same way to a surface (a forcein the z plane that does not affect migration in the xy plane nor the drag forces that oppose thatmigration) and thus 2D sheet formation is maintained. Additionally, the cells can polarize and exerta protrusive force (Factive) in the direction of open space, AECM. This force captures the ability ofmigrating cells to send out protrusions that bind to the free surface and allows the cell body to propelin the direction of the protrusion. I will refer to this force as adhesion to the ECM, where I considerthe ECM to be the free space that mimics the area created in experiments when a wound is made inthe epithelium (not to be confused to the previously mentioned 2D surface adhesion). I assume thatcells within the sheet can experience this force.The initial size and location of the cells is the same in all simulations. Fig. 7.1 shows an exampleof the initial condition for the special case where 40% of cells are red and 60% of cells are green.All simulations are initiated with 400 spherical cells and the radius of each cell is chosen from arandom distribution between (5,9) µm.I simulate growth as described in chapter 3, which implements cell growth of each of the threeaxes of the ellipsoid, at a constant rate. In the simulations shown, the growth rate is 0.003 min−1.If the cell experiences a large pressure from the surrounding cells the growth rate decreases andeventually stops (see details in Supporting Information A).7.4 Modelling resultsTo explore parameter space, I conduct numerous simulations varying parameters over a large range.I then refine the parameter space based on the observed morphologies of the epithelium and repeat85Figure 7.1: All simulations are initiated with 400 spherical cells distributed as shown in thisfigure. Here, 40% of cells are red and 60% of cells are green. The colors represent dif-ferent cell phenotypes that can have distinct adhesive properties. I explore how changingthese adhesion properties alters the distribution of the two cell phenotypes.the parameter sweep now in a refined range. Here, I show a subset of the results that demonstratehow the parameters influence the cell-sheet behaviour. The adhesion strength parameters are withina biologically realistic regime although, at present, there are no experimental measurements thatinform their choice.7.4.1 Changes in cell-cell adhesion, cell-ECM adhesion and the axis of cell divisionalter the leading-edge morphology of a cell-sheetI first study the leading edge morphology when all the cells have the same adhesive properties, seeFig. 7.2. The axis of cell division is chosen perpendicular to the net force acting on the cell. Whenthe cell-cell adhesion is strong enough the leading edge has a finger-like morphology, as seen inother models [83]. The onset of the instability that leads to the fingers emerging can also be causedby cell-ECM adhesion (recall that I define ECM to be the white open space in my simulations).Furthermore, my simulation results suggest that the cell-cell adhesion, cell-ECM adhesion and cellproliferation, strongly affect the leading edge morphology.When the cells have a strong adhesion to the ECM (AECM), relative to the cell-cell adhesion,fingers emerge, see Fig. 7.2 B. The results are shown at T=7 hours, which is characteristic of theleading-edge morphology throughout the simulation. Similarly, when the cell-cell adhesion is in-creased and the cells do not adhere to the ECM, fingering morphology emerges. Although bothadhesion properties lead to fingering morphology, the resulting shape of the fingers seem slightlydifferent. Most notably, the fingers that form when the cell-cell adhesion is strong are wider thanthose that form for strong cell-ECM adhesion. For strong cell-ECM adhesion, a leader-like pheno-type seems to emerge. The cells remain connected but this resamples the previously characterizedleader (or quasi-mesenchymal) phenotype [10].86A)Cell-Cell adhesionB)Cell-ECM adhesionFigure 7.2: Simulations of a wound-healing assay with a single cell type on a domain of size550 µm × 600 µm . Cells are moving towards the top of the figure. The results showhow A) increasing the strength of the cell-cell adhesion (0.1, 5.0 and 10.0 nN) or B)changing the strength of the cell-ECM adhesion (0, 1 and 1.7 nN), alters the leadingedge morphology. In A) there is no adhesion of cells to the ECM and the migration issolely driven by cell division. Results are shown at T = 7 hours.7.4.2 Initial ratios of cell types affects the rate of motion of the leading edgeI next focus on the emergence of the leader-cell phenotype by exploring how cells with differentadhesion properties organize and migrate in a cell-sheet. I conduct simulations with a mixture oftwo cell types, red and green cells, with different adhesion properties. I find that the adhesionbetween red cells, Arr, between green cells, Agg, and between red and green cells, Arg, all influ-87ence the distribution of the two cell types within the sheet. This result is reminiscent of predictionfrom Steinberg’s experiments of differential adhesion based cell sorting within aggregates [87]. Forexample, when the Arg is big relative to Arr and Agg a pattern that resamples a checkerboard is ob-served, see Fig. 7.3. This is to be expected since the adhesion strength is such that the cells preferto be surrounded by cells of the other color.Figure 7.3: Results from a simulation as in Fig. 7.2 but starting with 40% red cells and Arr =Agg < Arg. As is to be expected, the cells prefer to be in contact with cells of the oppositecolor resulting in a pattern that resembles a checkerboard. Results are shown at T = 5 hrswith Arr = Agg = 0.25 nN and Arg = 2.5 nN.In Fig. 7.4 the adhesion between cells is such that Arr > Agg > Arg and only the red cells adhereto the ECM, AECM. In these simulations, the red cells take over the leading edge of the sheet. InFig. 7.5 the position of the leading edge of the sheet is tracked in each of the simulations fromFig. 7.4. As the proportion of red cells is increased the cell-sheet seems to move more slowly. Atfirst, this result might seems surprising because the red cells adhere to the ECM. However, the redcells also have a stronger adhesion, meaning that they exert a greater pressure on their neighbourswhich slows down the cell growth, as explained in section 7.3.7.4.3 Cell-cell adhesions affect relative distributions of cell types and leading edgemorphologyFinally, I investigate how changes in adhesion properties for these two cell types alters their distribu-tion. Fig. 7.6 shows simulation results where the adhesion between red and green cells is increasedfrom left to right and the adhesion between red cells and the ECM increases from top to bottom. Ifthe AECM is strong enough, the red cells start breaking away. The Arg seems to have a lesser effecton the leading edge morphology, although the middle row shows that when Arg is large the red cellsdo not break away.88Initial fraction of red cellsFigure 7.4: Simulation setup as in Fig. 7.2 but the initial fraction of red cells, from left toright, is 0%, 20%, 40%, 60%, 80% and 100%. Only the red cells adhere to the ECMand Arr > Agg > Arg. The red cells take over the leading edge and as the proportion ofred cells increases, the sheet migrates more slowly. Results are shown at T = 7 hours.Parameters: Arr = 1.0 nN, Agg = 0.6 nN, Arg = 0.2 nN and AECM = 1.5 nN.Figure 7.5: Simulation results for Arr > Agg > Arg showing the location of the leading edge ofthe wound in the simulations with different initial ratio of red cells from Fig. 7.4. Thewound closure seems to slow down as the initial ratio of red cells is increased.7.5 Wound healing assay experimentsBased on the simulation results, the Roskelley lab has been conducting wound healing experimentswith cells that have different expression of cadherin adhesion molecules. The experiments are doneusing EpRas cells, which are pre-malignant mammary cell that have epithelial characteristics. Thefollowing three EpRas phenotypes expressing different adhesion molecules are used:89A ECMArgFigure 7.6: Simulation setup as in Fig. 7.2 but with 60% initial ratio of red cells. Each panelshows how the different adhesion between red and green cells, Arg, and red cells andECM, AECM, changes the leading-edge morphology. The results show that increasingArg makes it harder for the red cells to reach the leading edge whereas increasing AECMenables the red cells to take over the leading edge and eventually break away from thedensely packed sheet. Parameters: Arg = 0.0, 0.7, 1.0 nN and AECM = 0.0, 1.3, 2.0 nN.901. EpE0, cells that only express epithelial cadherin (E-cad),2. Ep0N, cells that only express neural cadherin (N-cad),3. EpEN, cells that express both E-cad and N-cad.These different cadherin adhesion molecules subtypes are observed to cluster together. However,less is known about the heterotypic binding affinity (E-cad binding to N-cad). To visualize thedifferent cell types in the experiments the cells were labelled with either green fluorescent protein(GFP) or red fluorescent proteins (mCherry).7.5.1 Experimental resultsExperiments were conducted with the three cell types described previously and a mixture of twodifferent cell types, see Fig. 7.7. Results are shown 20 hours after the scratch wound is made inthe monolayer. Based on previous studies of cadherins, and our hypothesis that homotypic adhesionbonds are stronger than heterotypic bonds, the adhesion of the EpE0 and Ep0N cells likely comparesin the following way: AE pE0−E pE0 >AE p0N−E p0N >AE pE0−E p0N [72]. However, less is known aboutadhesion properties of a cell expressing both cadherin types, the EpEN cell phenotype. This suggeststhat there is a knowledge gap that the modelling can address.The experiment in the first panel of Fig. 7.7, is conducted with only EpEN cells. The invasionfront has small finger-like protrusions. In the second panel, only EpE0 cells are used in the ex-periment and the results look similar to the the first panel. The third panel shows experiment withonly Ep0N cells where the cells at the leading edge have broken away from the sheet and migrateas individual cells. This was unsurprising since cells that express N-cad tend to have mesenchymalcharacteristics. In the fourth panel, 60% of the cells are EpEN and 40% are EpE0. A couple ofsmall fingers emerge at the leading edge and the wound seems to close faster. Notably, the invasionfront was mostly EpE0 cells. In the last panel, 60% of the cells are EpEN and 40% are Ep0N. TheEp0N cells take over the invasion front and they migrate as individual cells. Behind the invasiveEp0N cells, long fingers of EpEN cells can be seen.7.6 Comparing wound healing simulations to experimental dataI initially hypothesized that the EpEN cells would have the lowest cell-cell adhesion strength andthat, consequently they would take on a leader-cell phenotype. This assumption was based on thetheory that the heterotypic adhesion strengths would not be as strong as the homotypic adhesion.Comparing the simulations to the experiments suggests that the opposite is true. In the experiments,the EpE0 and the Ep0N cells take over the leading edge.To compare the experiments to the simulations, the adhesion properties of the three phenotypesneed to be identified. The low red-ECM adhesion, AECM, should be attributed to a more epithelial91Only EpEN OnlyEpE0 OnlyEp0N60%EpEN (green) + 40%EpE0 (red)60% EpEN (red)+ 40%Ep0N (green)Figure 7.7: Experiments, conducted by Tak Poon and Alannah Wilson in the Roskelley lab,using epithelial cells with different expression of N-cad and E-cad. First figure: onlycells expressing both N-cad and E-cad, EpEN cells. Second figure: only cells expressingE-cad, EpE0 cells. Third figure: only cells expressing N-cad, Ep0N cells. Fourth figure:60% EpEN cells and 40% EpE0 cells. Last figure: 60% EpEN cells and 40% Ep0Ncells. Note that the EpEN cells are green in the fourth panel but red in the last panel.The results show that the leading edge morphology differs between these experiments.The bottom shows the enlarged view of the leading-edge morphology for the blue boxin the second panel and the purple box in the fourth panel. A typical cell diameter isapproximately 20µm.cell-line (cells expressing E-cad) whereas strong AECM should be attributed to mesenchymal cell-line (cells expressing N-cad). Fig. 7.8 shows that when AECM is small, the red cells do not take overthe leading edge. For an intermediate value of AECM, the red cells take over the leading edge but thecell-sheet remains intact. This result is comparable to the fourth panel in Fig. 7.7. For strong AECMthe leader-cell phenotype breaks away from the sheet, reminiscent of the experimental results witha mixture of EpEN and Ep0N cells, last panel in Fig. 7.7.The second row in Fig. 7.8 demonstrates that as AECM is increased the average displacementof red cells that start close to the leading edge is significantly higher than that of green cells. Thisis to be expected and is seen in the top row of the same figure. The displacement data also showsthat cells that start further from the leading edge, in a densely packed sheet, do not migrate far.The bottom row illustrates the cell-tracks for every seventh cell of the initial 400 cells. In future92AECMFigure 7.8: Simulations with varying adhesion strength between red cells and the ECM, AECM(as in Fig. 7.6). Second row: the total displacement of the initial 400 cells (in units of 10µm) with respect to their initial position along the y-axis (also in units of 10 µm). Thebottom row shows the tracks for every seventh cell out of the initial 400 cells. As AECMincreases the red cells take over the leading edge and eventually break free. The AECMvalues used in this figure were 0, 1.3 and 2.0 nN.studies, this type of track data can be quantified along with the velocities of the cells, to compare toexperimental data and identify differences between phenotypes.7.6.1 Cell tracking algorithmBefore my involvement in the project, cell tracking was a laborious manual process in the Roskelleylab. Together with a fellow student, Dhananjay Bhaskar, I designed a pipeline, using the opensource program CellProfiler, that automatically detects and track cells from microscopy images [8].93Detecting cell boundaries in such a densely packed environment is one of the most challengingproblems in image segmentation. Some ongoing issues in cell tracking include: a) identifyingdaughter cells after cell division occurs, b) the fluorescence intensity changing over the course ofthe experiment, and c) cells showing up in both the red and the green fluoresce channels becauseof their dense packing. The pipeline as well as the experiments, are being refined for continuoususe in the Roskelley lab. This involves changing the experimental setup to include a nuclear stain.The nuclear stain is stable and it would not change over the course of the simulation. The pipelinecould then be improved by using the nucleus to identify cells instead of using the cell boundaries asit currently does. This might also improve the accuracy of the calculation of cell area and perimeter,but those details will have to be analyzed when the data has been generated.Preliminary results from the cell tracking algorithmFig. 7.9 shows preliminary results from the cell tracking algorithm used on imaging data from theexperiment with 60% EpEN and 40% EpE0 cells in panels A)-C). A corresponding simulationwith Arr = Agg > Arg is shown in panels D)-F). These parameters showed a similar morphologyof the leading-edge of the epithelium as the experiment. In the experiments the cells of the samephenotype seem to cluster in the bulk of the epithelium. Simulations showed a similar results forArr > Agg > Arg and therefore I show the perviously mentioned case here to focus on the adhesionproperties between different cell phenotypes. The experimental images were taken every 5 min.Panels A) and D) in Fig. 7.9 show the shape of the leading edge in the experiment and simula-tion, respectively. Panel B) shows the number of EpEN (green dashed line), EpE0 (red dotted line)and total (blue solid line) cells as a function of time. The cell number seems to increase linearlywith time. The same trend is seen in the simulation results in panel D), although the slopes arenot the same. This will be adjusted in future simulations. However, the tracking algorithm and theexperiments are still being refined. Quantitative data of cell division will be used to set the divisionrate in the simulations.Fig. 7.9 C) shows the 10th (red dotted line), 50th (blue solid line) and 90th (green dashed line)percentile for the cell areas obtained from the tracking algorithm. Area is calculated as the numberof pixels that the cell occupies which can be converted into µm2. Interestingly, the 90th percentilewas substantially higher than the 50th percentile and increased over time. This may suggest thatthere is a phenotype with a larger cell area emerging. As previously mentioned, migratory cells tendto have long protrusions that may result in a larger cell area. Suggestive of the leader-cell phenotype.The corresponding simulation result in panel F) does not show this trend. In this simulation, thecells are represented by ellipsoids and thus the protrusive cell shape can not be captured with thismodelling framework.94A) D)B) E)C) F)Figure 7.9: Preliminary data from cell tracking algorithm for the experiment with 60% EpEN(green) and 40% EpE0 (red) cells, A)-C). The corresponding simulation with Arr =Agg >Arg is shown in D)-F). A) and D) show the shape of the leading edge from the experimentsand simulations, respectively. B) and E) show the number of cells as a function of time.Each Frame Number from the experiment corresponds to 5 min. Lastly, the 10th (reddotted line), 50th (blue solid line) and 90th (green dashed line) percentiles for C) thecell area from the experiments and F) the cell volume from the simulations is shown.Parameters: Arr = Agg = 2.5 nN, Arg = 1.25 nN and AECM = 1.5 nN.7.7 DiscussionIn summary, I conducted simulations to study the role of adhesion on the emergence of a leader-cell phenotype at the migration front of a cell-sheet. These simulations then guided the design of95experiments in the Roskelley lab. Together, the simulations and experiments will help elucidate therole of adhesion on the cell-cell and cell-ECM interactions that facilitate migration of a cell-sheetand the emergence of the leader-cell phenotype.I found that by changing cell-cell or cell-ECM adhesion strengths in the simulations, the shapeof the leading edge can go from flat to the emergence of finger-like protrusions and eventually tocells breaking away from the sheet, see Fig. 7.2 (research question 1). I also found that the axis ofcell-division plays an important role in the finger-like morphology (research question 2).There are some similarities between epithelial sheet finger formation and the physics of inter-faces between two flowing fluids. In Hele-Shaw fluids it has been shown that fingering instabilitiesoccur when the interfacial tension (T ) is large enough relative to the product of the finger velocity(U) times the displaced fluid viscosity (µ) (i.e. the instability occurs above a certain value of Uµ/T )[81]. The instability in my simulations is caused by changes in the adhesion properties of cells. Thestrength of this adhesion can be interpreted as changing the interfacial tension and the viscosity ofthe sheet. For large enough values of the adhesion there is an onset of fingering morphology of theleading edge. Hence, my results are seemingly analogous to the instability condition for fingers influids which would be interesting to study in more detail in future research. An alternate approachto studying fingering instabilities in a cell monolayer was proposed in [55], where the tissue edge ismodelled as an interface with curvature and surface tension. The interfacial shape is then related toforces of cellular motility, leading to distinct patterns.I simulated a cell-sheet consisting of two cell phenotypes with distinct adhesion properties. Itracked the location of the leading edge in simulations with various initial fractions of two cell types.I showed, that increasing the initial fraction of the population with higher adhesion strength resultedin a slower migration of the leading edge, see Fig. 7.5 (research question 3).Furthermore, I showed that having the two cell populations, with different adhesion properties,sufficed to preferentially get one cell-type to take over the leading edge of a migrating cell-sheet,see Fig. 7.6 (research question 4). Additionally, if the cell-ECM interaction was strong enough,the cells started breaking off. These simulation results suggest that the leader-like phenotype canoriginate from cells located inside a cell-sheet that have distinct adhesion properties that allow themto migrate to the leading edge of the cell-sheet.It should be clearly noted that signalling pathways and mechanical properties have not beeneliminated. These may include the environmental conditions at the leading edge facilitating phe-notypic changes that lead to leader cell emergence or mechanical forces from neighbouring cellsinfluencing a cell’s internal biochemistry and cell contractility. This work merely demonstrates theimportance of cell adhesion in the possible emergence of leader cells.The modelling framework guided the design of experiments being conducted in the Roskelleylab, that will advance this research. The Roskelley group has made three different EpRas phenotypesthat express either E-cad (EpE0), N-cad (Ep0N) or both (EpEN), see Fig. 7.7. Interestingly, the96EpEN phenotype that was presumed to have the lowest cell-cell adhesion and be the leader-cellphenotype, did not behave as expected. Specifically, it seems to have low adhesion to the ECM eventhough it expresses N-cad. When EpEN cells were mixed with EpE0 cells, the EpE0 cells seemto have more leader-like behaviour and they dominate the leading edge of the cell-sheet, see fourthpanel in Fig. 7.7 (research question 5). This result requires further experiments to test the adhesionstrength of the different cell types. The expression data for N-cad and E-cad does not reliably givethe cell-cell adhesion strength. Instead, an experimental setup that allows a direct measurement offorce, i.e. atomic force microscopy or cells migrating on a substrate where deformations can bemeasured to calculate mechanical forces, would be beneficial.Through the development of the cell tracking pipeline, quantitative data was extracted fromthe experiments, see preliminary results in Fig. 7.9. The Roskelley lab is working on improvingthe experimental setup to optimize the quality of the cell tracking. The quantitative data can helpcompare experimental and simulation results and estimate model parameters, such as proliferationrates. Additionally, Dhananjay Baskar has experimented with the use of algorithms, such as machinelearning and the K-means algorithm, to track and classify the distribution, shape and velocity ofthe different cells in the epithelium. The simulations showed that the leading-edge migration andmorphology as well as the relative location of the red and green cells within the sheet changed whentheir cell-cell adhesion properties were altered, see Figs. 7.3 and 7.6. The clustering algorithms mayprovide additional insights from the experimental data on the adhesion properties of the cells thatcan be compared to the simulation results in future work.Other future experiments that would help further elucidate the emergence of the leader-cellphenotype include:(i) experiments with various ratios of the two cell types,(ii) experiments with a mixture of the EpE0 and Ep0N cell types, and(iii) experiments that can directly measure the strength of the adhesion bonds.Here, I have presented a framework that combines experiments and simulations to study migra-tion of an epithelium and the adhesion characteristics of leader-cell phenotypes. The simulationsinformed the experimental design. The experiments were then used to infer adhesion strengths ofcells expressing different cadherins, based on the predictions from the simulations. This frameworkwill be used to test hypothesis regarding cell adhesion and the leader-cell phenotype origin.This research is one piece in the puzzle that can bring us closer to understanding the biophysicalproperties that enable cancer cells to invade into surrounding tissues, blood vessels and lymphaticsystem.97Chapter 8Coupling biophysical and biochemicalprocesses to study cancer-cell spheroids8.1 IntroductionIn the previous chapter, I studied mixtures of cells with different adhesion properties. I showed howvariations in these adhesion properties can alter the morphology of the leading edge of a migratingcell-sheet. In practice, how did the different cell phenotypes emerge? In this chapter, I will modela biochemical network that couples to the biophysics of the cell, to account for distinct phenotypicbehaviour. I will study cell-cluster morphology in a 3D experimental setup with a specific focus onhow migration and proliferation influence the morphologies. Elongated cluster morphologies aresuggestive of containing cells with invasive potential which plays a key role in cancer metastasis.I continue to collaborate with the Roskelley lab at UBC, where Marcia Graves and Erin Bellconduct experiments with 3D breast cancer cell spheroids. By embedding the spheroids in variousextracellular matrix (ECM) they observe drastically different spheroid morphologies for cell linesexpressing a specific protein called Podocalyxin (podo). Podo has been identified as an independentmarker for poor patient outcome and metastatic risk [86]. Furthermore, podo has been shown toexhibit abnormal levels in some breast, prostate, liver, pancreatic, and kidney cancers. It plays a keyrole in lumen formation which is important for normal function of kidneys and breasts. Other groupshave focused on the formation of lumens in small podo expressing clusters of kidney cells, MDCKcells [6, 103]. Here, I will focus more on the morphology of multi-cellular clusters of breast cancercells. A theoretical model is uniquely able to help experimentalists develop and test hypotheses anddesign experiments to elucidate how podo may contribute to the range of spheroid morphologiesobserved experimentally. The mechanical properties of the tumor micro-environment can changeover time and they can vary spatially [74]. Therefore, I will investigate how the invasive potentialof cancer cells is affected when the cell-ECM interactions change.98Section 8.2 lists the questions that guide this research. In section 8.3, I outline and analyze themodel that I will use to capture interactions between the cells and the ECM. Section 8.4 describeshow the HyDiCell3D will be used to couple the biochemistry and biophysics. Section 8.5 showssimulation results with various cluster morphologies depending on biochemical and biophysical pa-rameters. In section 8.6, I will compare the results from the simulations to the spheroid experimentsfrom the Roskelley lab. In section 8.7, I will suggest future directions for this research project.8.2 Research questionsOverarching questions that guide research in collective cancer cell migration include: How does theECM facilitate cell invasion and intravasation? How do cells within a cell-cluster communicate?How are biophysical and biochemical processes coupled? My research is guided by the followingkey questions:1. How do the interactions of cells with the ECM affect cell polarity and migration?2. How does the signalling pathway for podo, analyzed in individual cells, compare to multi-cellular simulations?3. What can be inferred about proliferation and cell-cell adhesion from cell-cluster morphology?8.3 Biochemical model representing the podocalyxin dynamicsHere, I propose a biochemical model for the internal signalling dynamics that influence the local-ization of podo. This biochemical model will later be linked to the biophysics of cells using HyDi-Cell3D. Specifically, I calculate the surface area of the simulated cells in contact with the ECM anduse it as a parameter in the biochemical model. Subsequently, the biochemical model determinesthe cell polarity axis, completing the interaction loop between the biochemistry and biophysics.8.3.1 A molecular switch model for podocalyxinThe Roskelley lab has shown that the localization of podo changes in the presence of a bead coatedwith fibronectin (FN), as shown in Fig. 8.1. FN is an ECM protein that can bind to integrins,trans-membrane receptors that facilitate cell-extracellular matrix adhesion. The figure shows thelocalization of podo (second panel) and β1-integrin (first panel). In the absence of the bead thepodo and β1-integrin are uniformly distributed throughout the cell, as seen in the first three panelsof the top row in Fig. 8.1. When the cell is brought into contact with the FN bead the podo localizesto the surface opposite the bead, second panel in the bottom row of Fig. 8.1, and the β1–integrinlocalizes to the surface in contact with the bead, first panel in the bottom row.Further research on the podo localization was carried out by Mostov et al. [6]. The authorsshowed how knocking out RhoA eliminated any abnormal localization of podo and reinforced lumen99Merge DICB1 Integrin Podocalyxin10um FN-Coated BeadB1 Integrin Podocalyxin Merge DICSuspensionIntegrin Engagement Initiates Podocalyxin SegregationMergeMerge (Z-Series)PodocalyxinB1 IntegrinAttached To FN 30minFigure 8.1: Snapshots from single cell experiments, using a FN coated bead, conducted inthe Roskelley lab. In the absence of the bead (top row) β1-integrin (shown in the firstcolumn) and podo (shown in the second column) are uniformly distributed throughoutthe cell membrane. However, when the cell comes into contact with the bead (bottomrow) the β1-integrin and podo localize to opposite ends. The last column shows themicroscopy images (Differential Interference Contrast, DIC). Figure courtesy of MarciaGraves.formation in MDCK cells (see illustration of normal localization of podo in Fig. 8.2). RhoA is anintracellular protein in the family of the small GTPases and plays an important role in cell polarity.It is known to be responsible for cell contractility and coupling of cell internal biochemistry andmechanical properties [78].In [6], a series of experimental manipulations were used to identify the signalling pathway thatinfluences the podo distribution within the cell. Based on this work, I extracted the simplified sig-nalling pathway shown in Fig. 8.3. When β1-integrins bind, downstream effectors including thefocal adhesion kinase (FAK) are activated causing inhibition of RhoA and thus contractility. Down-stream effectors of RhoA are responsible for inhibiting the endocytosis (which causes relocation)of podo. Thus, when multiple β1-integrins bind the endocytosis of podo is high, resulting in podobeing localized to the membrane that is not in contact with the ECM (normal polarity) allowinglumens to form. From this simplified signalling diagram, it is evident that RhoA plays an importantrole in the localization of podo.100Figure 8.2: Illustration of the podocalyxin relocation. In normal (apical-basal) cell polarity,the podo is located away from the ECM interface and the lumens can form. Invertedpolarity is when the podo localizes to the cell-ECM interface.8.3.2 RhoA dynamicsThe Keshet group has conducted multiple analytical and numerical studies of the dynamics of RhoAand shown how it affects cell polarity. This work includes well-mixed models of RhoA [42], inter-actions of RhoA with other polarity proteins including Rac and Cdc42 [24] and spatial models of thedistribution of the molecules within a single cell [54]. In this research, I adapt a simple well-mixedmodel for the dynamics of RhoA from [42], and use the RhoA concentration to infer the localizationof podo.In the classic well-mixed models for RhoA dynamics there is an active membrane-bound formof RhoA, A, and an inactive cytosolic form, I. The general form of this model is:At = f (A)I−h(E)A andIt =− f (A)I+h(E)A.(8.1)The two functions f (A) and h(E) are the rates of activation and inactivation, respectively, and theycapture the dynamics of switching between the states A and I.101RhoAROCKIFAK RacEzrin/NHERFPodoendocytosisoutside cellinside cellβ1-integrinFigure 8.3: Schematic of the podocalyxin signalling pathway extracted from [6]. podo andβ1-integrin are observed to localize to opposite areas (similar to the relationship betweenRac and Rho during cell polarization). The schematic shows that RhoA dynamics playan important role in the podo localization.A (Cell-ECM podo)I(Normal podo)E(ECM contact)bg훽훿Figure 8.4: Diagram of the simplified signalling pathway for the localization of podo. Thereis a transition between the A (podo localizing to the ECM interface) and I (the normallocalization of podo) states and a positive feedback for A. The ECM variable, E, willbe proportional to the surface area of a cell in contact with the ECM in multi-cellularsimulations.8.3.3 ODE model to describe podocalyxin localization within a cellHere, I build on the well-studied equations for RhoA from [42]. The model is outlined in Fig. 8.4and it captures the switching of podo between the cell-ECM interface, represented as A, and thenormal localization, on the membrane away from the ECM where the lumens would form, rep-resented as I. The positive feedback loop for A, f (A) is equivelent to the one used in [42] but Iadded the feedback from the ECM, in the inactivation rate h(E). This feedback is based on vari-ous experimental observations [6, 36, 103]. I used the following sigmoidal functions for the two102feedbacks:f (A) = b∗+β ∗An1+An,h(E) = g+δ ∗Em1+Em,(8.2)where b∗ and g are background activation and inactivation rates, respectively, β ∗ is the rate of thepositive feedback loop, E is the environmental parameter (in the multi-cellular simulations thiswill be proportional to the surface area of the cell in contact with the ECM) and δ ∗ is the ECM-amplified rate of podo relocation. In the biological system δ can be influenced by the availability ofβ1-integrin binding sites or other ECM factors that can alter the relocation of podo. The exponents,n and m, control the shape of the sigmoidal functions for the feedbacks. Note, that A and E havebeen normalized with respect to their inhibition level that yield half of the basal level (IC50), toobtain 1 in the denominator of the sigmoidal functions.In Eqn. (8.1) the total amount of podo is constant, AT = A+ I, and the transport time betweenthe states, A and I, is neglected. Thus, the system can be reduced to the following single ordinarydifferential equation (ODE):At = g[(b+βAn1+An)(AT −A)−(1+δEm1+Em)A], (8.3)where x = x∗/g for x ∈ {b, β , δ} are the unit-less rate parameters (b is the background activationrate, β is the positive feedback rate and δ is the ECM-amplified rate of podo relocation). Thisequation is known to have regions where bistability is observed [42]. Here, Eqn. (8.3) acts as abistable molecular switch between the podo state of the cell which will later determine the cellpolarity. Fig. 8.5 is a plot of the steady states of Eqn. (8.3) using δ as a bifurcation parameter. Thissystem has three different regimes, which I interpret as: I: podo localized to the cell-ECM interface,II: bistable region, III: podo localized away from the ECM interface.8.3.4 Sharp switch analysisThe powers in the two sigmoidal functions, n and m, determine the steepness of the feedback re-sponses, i.e. how quickly the feedback functions switch from low to high values. It is helpful toconsider the special case where (n,m)→∞, since the sigmoidal functions can be approximated withstep functions which greatly simplifies the analysis (similar to the analysis in section 6.4). In thislimit, the following piece-wise constant functions are obtained at steady state,1030 1 2 300.511.52I II IIIFigure 8.5: Bifurcation diagram showing the three regimes of behaviour using δ , the ECM-amplified rate of podo relocation, as a bifurcation parameter: I: podo localizes at the cell-ECM interface, II: bistable region, and III: podo localizes away from the ECM (lumenscan form). Here, h(E) = g+ δ and the parameters used are: b = 0.2, β = 4.0, n = 3,m = 2, AT = 2 nM, g = 0.5 s−1.F (A) =b+β , A 1,b, A 1,H (E) =g+δ , E 1,g, E 1.Eqn. (8.3) can now be written in the following form,At = g [F (A)(AT −A)−H (E)A] . (8.4)This equation has four different steady states, namely:A1 =(b+β )ATb+β+g+δ A 1, E 1,A2 =(b+β )ATb+β+g A 1, E 1,A3 = bATb+g+δ A 1, E 1,A4 = bATb+g A 1, E 1.(8.5)The A1 and A2 steady states represent high A (podo located away from neighbouring cells) and A3and A4 are low A. The regimes where these steady states exist are shown in Fig. 8.6, using δ and104AT as bifurcation parameters. The δ parameter gives insight into how the ECM influences the cellbehaviour and the total amount of podo, AT , varies between cell types which may be responsible forthe differences in spheroid morphologies. The steady states in the regions labelled in the figure are:(a) A3 and A4, (b) A2, A3, A4, (c) A1, A2, A3, A4, (d) A1, A2, A3, (e) A1 and A2. The bistability occursin the shaded regions (b) and (c) for E 1 and in shaded regions (c) and (d) for E 1.Figure 8.6: Two-parameter bifurcation plot for δ , the ECM-amplified rate of podo relocation,and AT , the total amount of podo. The figure shows the five different regimes obtainedfrom the sharp switch analysis in Eqn. (8.5). The lines from left to right are AT = 1+gb+β(red, dash-dotted line), AT = 1+g+δb+β (blue, solid line), AT = 1+gb (red, dashed line) andAT = 1+g+δb (blue, dotted line) . The steady states are: (a) A3 and A4, (b) A2, A3, A4, (c)A1, A2, A3, A4, (d) A1, A2, A3, (e) A1 and A2. Bistability occurs in the shaded regions ((b)and (c) for E 1 and (c) and (d) for E 1).This theoretical consideration of the localization of podo emphasizes the rich dynamics of themodel. In practice, the interactions with the environment, δ , can vary both in space and time.Consequently, there would be a transition between the possible steady states.The model can be related to the experiment with the FN bead in Fig. 8.1. In the absence ofthe bead, there are few available β1-integrin binding sites, captured with the δ parameter beingsmall. This causes A to be large and podo localizes to the the cell-ECM interface. As the FN bead isbrought into contact with the cell, the δ parameter increases and the podo localizes away from thebead, representing the ECM interface. These results show the dynamics of the model in single cells.I will now explore how cells in a multi-cellular system will behave when such internal signalling isincluded in their dynamics.1058.4 Multi-cellular simulations coupling biochemistry to biophysicsI ask how a multi-cellular system behaves when the polarity of each cell is determined based onthe model derived in the previous section. I use the HyDiCell3D framework which is capable ofcoupling biochemical and biophysical properties in a multi-cellular system. Specifically, I use theODE in Eqn. (8.3) to determine the polarity axis for each cell. The ECM variable, E, is proportionalto the amount of surface area that is not in contact with other cells. I also ask how observations ofthe multi-cellular system can be used to infer the properties of individual cells.Each simulated cell can be in one of the following two phenotypic states based on the internalbiochemistry:Green cells: when A > Ath the podo is assumed to be localized to the cell-ECM interface (awayfrom the neighbouring cells) and the polarity of the cell is inverted. The polarity axis is determinedto be pointing away from the neighbouring cells.Red cells: A≤ Ath, podo is assumed to localize towards neighbouring cells and the cell has normalpolarity. The polarity axis is chosen in the direction of the neighbouring cells.8.4.1 HyDiCell3D simulation outlineThe following is an outline of the simulations.• The simulations are initialized with 400 single cells distributed with a set distance, dspread ,between them. To obtain initial conditions comparable to the spheroid experiment, I first carryout the simulations with only the mechanics (no signalling) until clusters form (simulationtime is 8 hours). All the cells start with the inverted polarity meaning that the podo is locatedat the cell-ECM interface (green cells), see Fig. 8.7.• Next, I activate the dynamics from Eqn. (8.3) that captures the internal biochemistry in eachcell. Table 8.1 lists the baseline parameters used in the simulations.• The E parameter in Eqn. (8.3) is proportional to the surface area of the cell that is not incontact with other cells, and assumed to be interacting with the ECM. See further details inSupporting Information A.• In each time step, the value of A from Eqn. (8.3) that represents the localization of podo,determines the polarity of the cell, as explained earlier.• The polarity determines the direction of the cell-ECM adhesion force. The direction of thisforce can alter the contact area of cells with the ECM, the parameter E, that in turn influencesthe internal biochemistry. Thus, the interaction loop between the internal biochemistry andthe mechanics is completed, see illustration in Fig. 8.8.106The parameters in Table 8.1 for the biochemical model are not based on literature since thismodel is phenomenological. The cell-cell adhesion and cell-ECM adhesion strengths are chosen tobe within a biologically relevant parameter regime, which has been reported in the range of 0.1-20nN [7, 18, 51]. However, I do not have direct measurements of these forces for the specific celllines and ECM used in this research. There is some evidence that podo decreases the strength ofcell-cell adhesion bonds [90]. Therefore, the podo expressing cells may have adhesion strengths onthe lower end of this range.Pre-clustering t = 0 t = 90 min t = 30 hrsFigure 8.7: Left to right: time series from a simulation where the ODE from Eqn. (8.3) issolved in each cell to determine the cell polarity. Red cells: normal polarity. Greencells: podo localizes away from neighbouring cells. First panel: the initial spreading ofindividual cells in the simulations. Second panel: the clusters have formed and now theinternal biochemical signalling is activated. The last two panels show how the clustermorphology changes over time. Parameters as in Table 8.1.8.4.2 Cluster morphology quantificationIn order to quantitatively compare the results from the simulations to the experiments, I calculate anElongation Index for clusters, both from experiments and simulations. I use the Image ProcessingToolbox in Matlab for this calculation. I define the Elongation Index as the ratio of the major tominor axis of an ellipsoid that best matches the cluster. For a spherical cluster the elongation indexwould be close to 1 and as a cluster elongates, the index increases.This index measurement works well for elongated morphologies but for clusters that have ir-regular shapes, for instance with many protrusions, a different index calculation would be required.Furthermore, the spheroids used in the experiments are 3D, I only obtain 2D phase images of the re-sults. Therefore, I quantify a 2D property of a 3D object, restricting the quality of the measurement.For consistency, I also calculate the elongation index of 2D images of the 3D simulations.107Intracellular signallingCell polarityCluster morphologyAIECMFigure 8.8: Illustration of how the intra-cellular signalling influences the cell polarity whichleads to a change in the tissue morphology. Eventually the new tissue morphology altersthe variable E, cell surface area in contact with the ECM in Eqn. (8.3), that influencesthe internal biochemistry. This completes the indirect signalling loop between the bio-chemistry and biophysics in the simulated cells.8.5 Simulation resultsSimilar to the parametric exploration procedure I used and described in the previous chapter, Iperform a number of simulations with various parameters. The figures shown in this section capturethe breadth of results for both cluster morphologies and cell polarity within a cluster.8.5.1 The biochemical network influences cell polarityI conduct simulations, using the setup described in the previous section, and analyze the morphologyand distribution of cell types in the resulting clusters. I ask how the biochemical network needs tocouple to the cell mechanics to obtain the morphologies observed in experiments. I find that it issufficient for the state of the cell, A, to determine the cell polarity. More specifically, this assumptionleads to the formation of clusters that, in certain parameter regimes, are loose and appear to havelumens. For other parameters, the results show elongated cluster morphologies with a mixture ofcells with different polarities (red and green cells).Because of the anti-adhesive property of podo, I experimented with having the cell state, A,108Table 8.1: Summary of the baseline model parameters for simulations of the podocalyxinlocalization in multi-cellular clusters.Parameter Definition Units Valuedspread Initial spacing between cell centres µm 17Arr Adhesion between red cells nN 0.25Agg Adhesion between green cells nN 0.25Arg Adhesion between red and green cells nN 0.25AECM Adhesion of red cells to the ECM nN 0.2δ ECM amplified rate of podo relocation – 4β Rate of the positive feedback – 4g Rate of switching from A to I 1/s 0.5b Rate of switching from I to A – 0.14m The power in the positive feedback loop – 3n The power in the feedback from the ECM – 2AT Total amount of podo nM 2Ath Threshold value of podo for polarity nM 0.6·ATinfluence the cell-cell adhesion strength as well. I found that this added complexity did not resultin any new cluster phenotypes in the simulations. Therefore, I decided to not include it in this workbut this assumption can be revisited for example if the lumen formation is to be studied in greaterdetail.Lumen formation Lumen formation is an active process where either, the inner cells die to formthe lumen (as is the case for MCF-10A mammary cells) or they are repositioned due to intercellularcavitation (i.e. MDCK cells). The MCF-7 cells used in this research tend to form tumors andnormally they do not form lumens. The podo has anti-adhesive properties which can influencelumen formation. In the simulations, by using the biochemical model to only infer the polarity axisfor each cell, the lumens could form given the right balance between cell-cell adhesion and cell-ECM adhesion. As the cell-cell adhesion is increased the results go from loose clusters to cohesiveclusters with both normal and inverted polarity (red and green cells, respectively), see Fig. 8.9.Fig. 8.10 A) shows simulations with small cell-cell adhesion (left) that form loose clusters and largecell-cell adhesion (right) that form cohesive clusters. Fig. 8.10 B) demonstrates that given the rightset of parameters and small enough clusters, lumens are observed in the simulations. The first panelshows a top-down view (along the z-axis) of a cluster. The second panel shows a cross-section, inthe x-y plane, through the middle of that cluster. I show the same cross section again but with thex-y plane rotated 55◦ in panel three and 87◦ in the last panel.109Figure 8.9: Simulation results showing how varying the strength of the adhesion between cellsalters the cluster morphology and distribution of cells with different polarity (red andgreen cells). For low cell-cell adhesion (first panel) the clusters are loose and most ofthe cells have normal polarity (red). For high cell-cell adhesion (last panel), elongatedclusters with partial inverted polarity are observed. Agg = Arr = Arg = 0.25, 0.35 and 0.5nN, other parameters as in Table 8.1.8.5.2 Variations in the total amount of podocalyxin in each cell alter cell polarityThe analysis of the biochemical model, see Fig. 8.6, shows that variations in the total amount ofpodo, AT , change the possible steady states describing podo localization. To compare the resultsfor single cells to multi-cellular simulations, I vary the total amount of podo available in each cell,see Fig. 8.11. For low total amounts of podo the cells have normal polarity and loose clusters areobserved. For intermediate amounts of podo, in the bistability regime, cell-clusters are observedwith both normal and inverted polarity. For large total amounts of podo all the cells have invertedpolarity.To compare the trends of the bistability zones from the analysis in individual cells to the multi-cellular simulations, I repeated the simulations with various parameters δ , see Fig. 8.12. The sizeof the bistability zone expands as the parameter δ is increased. This trend is the same as predictedin Fig. 8.6.Next, I will assess how the polarity reversal leads to invasion of cells into the ECM (character-ized by changes in cluster morphology).8.5.3 The micro-environment influences the cluster morphologyBefore I investigate the morphology of clusters based on the localization of podo, I obtain the mor-phology of clusters from the baseline simulation where the cells do not express any podo. Turningoff podo in the simulations means that I simulate only the mechanics of the cell interactions without110A)B)Cluster Cross section (CS)CS rotate 55° CS rotate 87°Figure 8.10: A) simulations with different density of cells within a cluster. Left panel: cell-celladhesion is 0.25 nN and AT =2 nN, the cluster is loose. Right panel: cell-cell adhesionis 0.5 nN and AT =1 nN, cell-ECM resulting in a cohesive cluster. Other parameters asin Table. 8.1. B) Simulation results showing a slice through a loose cell-cluster (firstpanel), from three view points, to demonstrate that in certain cases, when the cluster issmall enough, there is a lumen forming in the center. Second panel is a cross sectionthrough the middle of the cluster, in the x-y plane. The same figure is shown with thex-y plane rotated 55◦, third panel, and 87◦, last panel.any internal biochemistry. The results in Fig. 8.13 show simulations without podo (grey cells) onthe left and with podo (red and green cells) on the right. The control cell-clusters are cohesive andhave a spherical morphology whereas the podo expressing cells have a mixture of cell polarities andit certain parameter regimes the cluster morphologies are elongated. In the next section, I calculatethe elongation index for the different simulations to quantitatively compare the simulations.I next investigate how varying the ECM input changes the results in my model. There aretwo important parameters in the model for the ECM interactions: δ , the rate associated with theECM input on the biochemical network and AECM the strength of the adhesion between red cellsand the ECM. Fig. 8.14 shows how varying these two parameters altered the morphology ofthe cell-clusters. In the cluster morphology feature space, I classify the following three featuretypes: i) normal polarity, ii) partially inverted polarity and iii) inverted polarity. These three clusterphenotypes are related to the three steady state regimes for the biochemical model in Eqn. (8.3).As δ increases, from left to right, or AECM increases, from top to bottom, the number of cellswith normal polarity, red cells, increases and the loose clusters are observed. The normal polarity111AT=1 AT=2 AT=20Figure 8.11: Simulation results showing how varying the total amount of podo, AT , changesthe morphology of the clusters. When AT is increased the results switched from normalpolarity, to bistability, where both steady states are present, and eventually to invertedpolarity with podo at the cell-ECM interface. This is comparable to the analysis of thebiochemistry in individual cells illustrated in Fig. 8.6. Other parameters as in Table 8.1.phenotype is observed when both δ and AECM are large and the inverted polarity is observed whenδ is small. For large values of AECM the clusters start to disassociate and individual cells break awayfrom the clusters.8.5.4 Axis of cell division plays a key role in cluster morphologyThe doubling time of the podo expressing breast cancer cells has not been measured in the spheroidexperimental setup but it is estimated to be of the order of 12-24 hours. Since the spheroid exper-iments are conducted over a few days, it is safe to assume cell division occurs. I next ask whateffect it will have in the simulations to include cell division. I find that implementing cell prolif-eration in the simulations naturally results in bigger clusters and the density of clusters increases,see Fig. 8.15. Additionally, more of the outer-most layer of cells in the cluster have normal polarity(red cells).Interestingly, the axis of cell division greatly alters the simulation results. The middle row inFig. 8.15 shows the results when the axis of cell division is chosen perpendicular to the polarityaxis of the cells. The bottom row in the same figure shows the results when the cell division axis ischosen at random from the three primary axes of the ellipsoids. The results demonstrate that boththe morphology and the cell polarity vary between these two simulations. This result suggests thatthe axis of cell division may play an important role in the elongation of the spheroids and that mightbe due to the matrix interactions altering the polarity of the cells.112Figure 8.12: The size of the zone where bistability occurs changes for various values of thetotal amount of podo, AT and the strength of the environmental feedback, δ , as theanalytical results for individual cells in Fig. 8.6 predicted. This plot shows the fractionof cells with normal polarity (red cells) as a function of AT for four different values ofδ . The fraction of cells with normal polarity can be used to infer the clusters that havenormal polarity. Below the brown horizontal line, only 95% of the cells have invertedpolarity and above the purple horizontal line 95% of cells have normal polarity. Thezone between these horizontal lines is where partially inverted polarity can occur. AsI predicted from the analysis, by increasing δ the bistability zone (partially invertedclusters) expands. Other parameters used are as in Table 8.1.8.6 Comparing multi-cellular simulations with podocalyxin model toexperimental resultsI compare my simulation results to experimental data from the Roskelley lab. The experimentsare conducted with either control MCF7 breast cancer cells or MCF7 cells with up-regulated podoexpression, MCF7-podo [37]. First the cells are grown overnight on a plate so that 3D spheroidsform. Next, the spheroids are overlaid with either Matrigel or collagen-I, see results in Fig. 8.16.Little difference is observed between the morphology of the clusters from these two cell types whenoverlaid with the Matrigel. However, when they are overlaid with the collagen-I, the MCF7-podospheroids invade into the collagen, resulting in elongated cluster-morphologies. Meanwhile, theMCF7 control cells in collagen-I remain spherical.113Control cells Cells interacting with ECMFigure 8.13: Simulation results showing how the cluster morphology changes when the bio-chemical network is activated. The control cells, left panel, do not have the biochemicalnetwork and the spheroids have a cohesive rounded morphology. When the signallingnetwork is activated, right panel, the clusters consist of cells with different polaritywhich can alter the cluster morphology. The cluster shapes will be quantified later, seeFig. 8.18. The parameters in the right panel are as in Table 8.1.By comparing the experimental results to the simulations in Fig. 8.14, the cell-ECM interactionscan be inferred. The δ parameter, the ECM amplified rate of podo relocation, can be assumed tobe higher in Matrigel than in collagen-I. Matrigel has more laminin which has been shown to playan important role in the redistribution of podo [103]. Thus, the MCF7-podo cell line overlaid withcollagen-I can be compared to the panels in the first column of Fig. 8.14. The results are in goodagreement with the observed elongated morphology. Similarly, the MCF7-podo cell-line overlaidwith Matrigel can be compared to the bottom right section of the same figure (both δ and AECMlarge) where the spherical morphology and normal polarity are observed. The control MCF7 cellshave a cohesive spherical morphology reminiscent of the baseline simulation in the left panel ofFig. 8.13.8.6.1 Quantitatively comparing cluster morphologiesIn order to quantitatively compare the results from the simulations to the experiments, I calculate theelongation index for both the experimental and simulation results. First, I segment the image usingthe image processing toolbox in Matlab to identify cluster boundaries and eliminate clusters that are114휹A ECMFigure 8.14: Simulation results showing how varying the parameter δ and the strength of theadhesion force to the ECM, AECM, alteres the morphology and size of the clusters.Results are shown at t=3 days. Parameters: δ = 1.0, 2.0, 10.0, 20.0, AECM = 0.0, 0.2,0.4, 0.5 nN and other parameters as in Table. 8.1.too small as well as clusters that touch the boarder of the image. Fig. 8.17 shows the phase imagesfrom the experiments and the corresponding segmented image. The elongation index calculatedfrom the simulations follows the same trend as the one calculated for the experiments, Fig. 8.18.Interestingly, the simulations indicate that the addition of proliferation, with the axis of cell divisionchosen perpendicular to the polarity axis, does not significantly alter the elongated morphology of115t = 0 t = 1 day t = 3 daysFigure 8.15: Varying cell proliferation alters the cluster morphologies. Top row: no prolifer-ation. Middle row: the axis of cell division is chosen perpendicular to the polarity axisof the cell. Bottom row: the axis of cell division is chosen randomly from the threeprimary axes of the simulated ellipsoid. The cell doubling time is 23 hours. The resultsshow that proliferation changes the density and size of the clusters, as expected, as wellas the polarity of the cells. δ=1.0 min−1, AECM=0.1 nN and other parameters as inTable 8.1.the clusters, see purple circles in Fig. 8.18. However, qualitatively, it is evident in Fig. 8.15 that theproliferation alters the polarity of cells within the cluster.116Figure 8.16: Phase image results from experiments using the control MCF7 breast cancer cell-line, left column, and a MCF7-podo breast cancer cell-line with up-regulated expressionof podo, right column. The top row shows the spheroids before they are overlaid withECM. The clusters that are overlaid with Matrigel are shown in the middle row andthose overlaid with collagen-I, are in the bottom row. Little difference is observed be-tween the two cell types in Matrigel whereas in the collagen-I the MCF7-podo cell-linehas a significantly elongated morphology compared to the control cells, as is quantifiedin Fig. 8.18. Figure courtesy of Dr. Marcia Graves.8.6.2 Comparing different extracellular matricesThe two parameters that are altered in Fig. 8.14 have a physical interpretation. As stated before,the δ parameter captures the ECM induced redistribution of podo which can be caused by lamininor β1-integrin binding sites. The availability of laminin is higher in Matrigel than in collagen-I. Additionally, β1-integrin blockers were used on MDCK cells in [6] which changed the clustermorphology from normal polarity to partially inverted polarity. This is equivalent to the results in117Figure 8.17: The phase images, first column, and the corresponding segmented image, secondcolumn. The segmentation is conducted in Matlab and all clusters that touch the pictureboarder or that are too small to be multi-cellular clusters, are removed.Fig. 8.14 where moving from right to left in the bottom two rows changes the clusters from normalpolarity to partially inverted polarity.8.6.3 Force generation of cells in contact with the extracellular matrixThe cell-ECM adhesion, AECM, which increases from top to bottom in Fig. 8.14, depends on internalforce generation, connected to the actomyosin network. Experiments have been conducted to blockthe actomyosin contractility, using blebbistatin, in MCF7-podo cells overlaid with collagen-I [37].The elongated cluster morphology disappeared in these experiments. This should be compared tothe simulation results in Fig. 8.14 going from the bottom to the top. Since the cells were embeddedin collagen-I, the δ parameter is small. These simulations do not seem to agree with the experi-mental results. However, some research has shown that myosin contractility also plays an importantrole in epithelial cell adhesion [104]. Furthermore, myosin can play a role in the localization ofβ1-integrin and cell polarization. Fig. 8.9 shows that decreasing the cell-cell adhesion parametercan change the resulting cluster morphology from elongated to spherical. I thus conclude that theblebbistatin experiment can not be compared to simulations where the cell–ECM interactions de-118A) B)ControlMatrigelCollagen11.522.533.5ControlMatrigelCollagenProliferation11.522.533.5Figure 8.18: Using the image processing toolbox in Matlab, I calculate the ratio between themajor and minor axes of an ellipsoid that best matches the clusters from the experimen-tal and simulation data. A) the elongation indices for the three experimental setups,namely the control cells in collagen-I and the MCF7-podo cells cultured both in Ma-trigel and collagen-I. B) the elongation indices from the simulation results. The prolifer-ation simulations are conducted with the same parameters as the collagen-I simulation.crease and that the myosin network plays a complex role in the cluster dynamics which requiresfurther analysis.8.6.4 Manipulating the RhoA pathwayLastly, Fig. 8.11 can be compared to experiments where the RhoA has been knocked down, as wasdone in MDCK cells in [6]. The experiments showed that RhoA knock-down abolished the par-tially inverted cluster polarity. Recall that I use a dynamical model for RhoA (outlined in Fig. 8.4)to infer the location of podo. I now interpret the results in Fig. 8.11 for AT being the total amountof RhoA. These simulations show that when RhoA decreases from 2 to 1nM the results are equiv-alent to the experimental observations with partially inverted polarity. Additionally, the previouslymentioned blebbistatin experiments block actomyosin contractility (presumably mediated throughRhoA) which abolished the elongated morphologies [37]. Those results align well with the decreasein total RhoA in Fig. 8.11.1198.7 Future work8.7.1 Modelling alignment and density of the ECMA future improvement of this cancer-cell spheroid study should involve investigating the structureof the ECM and how proteins that the cells secrete to alter their micro-environment (matrix metallo-proteinases, MMPs) influence collective cell migration. Certain MMPs have been observed to playimportant roles in cell invasion and metastasis [62]. This research could be done with simulationsof cell clusters that influence the ECM density, the collagen fibre orientation, or both. Fig. 8.19demonstrates an example of such interactions in the HyDiCell3D where the environment (whichhas been discretized) starts with random alignment and density that the cells can alter during migra-tion. In this example, the red cells degrade the ECM at a constant rate, they can also release newECM and they alter the alignment in each discretized voxel by averaging the directional vector ofthe ECM with the cell’s direction of migration.Figure 8.19: Simulations using a field to represent the density and random orientation of theECM (collagen fibers). These cells do not have the biochemical model for the podolocalization. They have cell-cell and cell-ECM adhesion and they polarize based on thefield orientation. From left to right: 0, 1, 2 and 3 hrs. The simulations are conductedwith 57 green cells and 7 red cells. The red cells interact with the ECM by altering thedensity and fibre orientation.8.7.2 Cluster migrationIn experimental work by Bryant et al. [6], the migration of podo expressing clusters in differentconditions was tracked. The authors observed that clusters with partially inverted polarity of podohad greater migration than clusters with normal polarity. Fig. 8.20 shows the paths of every othercell in two different simulations, one where clusters have partially inverted polarity, top row, andone where the clusters have normal polarity, bottom row. The cluster migration is governed bythe cell-ECM adhesion and the cell polarity. Blue paths represent slow movement and red pathsrepresent fast movements. These results show that the clusters with normal polarity do not move120far, bottom row, contrary to the partially inverted clusters which move quite a bit, top row. Futureresearch should include quantifying and analyzing this migration both in the simulations and theexperiments in the Roskelley lab.Figure 8.20: Simulations with partially inverted polarity, top row, and normal polarity, bottomrow. The first column shows a snapshot at the end of the simulation and the rightcolumn shows the tracks of 200 out of the 400 simulated cells. Blue tracks representslow movement and red tracks represent fast movement.8.7.3 Impeding the RhoA activation experimentallyThe Roskelley lab is experimenting with using cell mutants with podo that cannot interact withEzrin/NHERF, a downstream effector of RhoA (see Fig. 8.3). The endogenous wildtype podo hasbeen knocked-out in these mutated cells. The idea is to use these mutants to elucidate the role ofRhoA during morphological changes of the cell clusters.1218.8 DiscussionI studied the morphological changes of 3D cell spheroids. I focused on the role of Podocalyxin(podo) and investigated spheroid morphologies and polarity of the cells within a spheroid. Using asimplified model motivated by the dynamics of RhoA, I determined the polarity of individual cells.The signalling network acts as a bistable molecular switch and I determined how the influence of theECM can act as a bifurcation parameter between the three states; normal polarity of podo, invertedpolarity of podo and a bistable state (research question 1). I implemented this signalling modelin the HyDiCell3D framework and analyzed the spheroid morphologies. I found good qualitativeagreement between the analytical predictions for the bistability region, Fig. 8.6, and the multi-cellular simulations Fig. 8.12 (research question 2).I used experiments from the Roskelley lab with various cell-lines and ECM conditions to com-pare to my model. Specifically, I interpreted the various cell interactions with the ECM as differentadhesive strengths in my model or parameters in the biochemical network. The model is able tocapture the differences in experiments with spheroids embedded in either Matrigel or collagen-I, asis quantified using an elongation index to describe the morphology of the spheroids, see Fig. 8.18.Fig. 8.21 summarizes the different morphological phenotypes observed in experiments and cap-tured in the simulations. Interestingly, I found that including cell proliferation in the simulationsdid not have a significant effect on the elongation index when the axis of cell division was cho-sen perpendicular to the cell polarity axis. However, when the axis of cell division was chosenrandomly the resulting cluster morphologies, as well as the polarity of the outer most cells in theclusters, changed (research question 3). This model prediction would be interesting to explore inexperiments especially since abnormal cell division plays a role in the onset of a number of diseases,including liver disease and cancer [67].Following is a list of experimental and theoretical work that I suggest as future directions andvalidations of this research:I Repeating experiments with various densities of collagen-I or different amounts of availableβ1-integrin binding sites to study the range of cluster morphologies. Additionally, the den-sity of collagen-I could be kept constant while the matrix stiffness is altered (by combiningcollagen-I with cross-linkable polyacrylamide containing hydrogels). Experimental observa-tions have been made that suggest that lumens are unable to form in stiff collagen matrices[74]. My simulation results suggest that a range of spheroid phenotypes could be obtained byaltering these environmental conditions, see Fig. 8.14.I Repeating the experiments with lower density of initial clusters. In previously published workfrom the Roskelley lab using 4T1 cells, some of the clusters merged over time, contributingto the morphological elongation of the clusters [37]. These experiments could elucidate howmuch of the cluster elongation is attributed to cluster merging versus cluster migration.122 147 147 Podo / NucleiSimulationsExperimentsFigure 8.21: Comparing the cluster morphologies observed in experiments to simulation re-sults. The three phenotypes are: cohesive clusters (left, MCF7 control), elongated clus-ters (middle, MCF7 podo) and clusters with lumen (right, EpH4 cell-line in Matrigel).Experimental figures courtesy of Dr. Marcia Graves.I Conducting simulations where the axis of cell division could be detected, for example throughvisualizing the microtubules. My simulations suggest that the morphologies of the spheroidsare greatly affected by this division axis, see Fig. 8.15.I Exploring the interplay between the podo and actomyosin networks and how that may influ-ence the spheroid morphologies. One possibility is to conduct experiments with mutant cellsthat do not bind Ezrin/NHERF (section 8.7.3).I Developing a spatial experimental and computational model for the localization of podo andβ1-integrin.I Repeating the simulations with a number of random initial conditions to extract statisticalinformation.I Replacing the elongation index with a shape classification that can capture more complexmorphologies. This may include using supervised machine learning algorithms or a differentshape measurement that could capture size and shape of protrusions extended out from thespheroid.123I Investigating the structure of the ECM and the migration of the spheroids as described insection 8.7Cells are known to respond to physical changes in the ECM through both internal signallingpathways and mechanical forces. This research highlights that the coupling between mechanicsand biochemistry is an important research area that can provide novel insights into the invasion ofcancer cells. Here, I developed a framework that can be used to test hypotheses regarding mechano-chemical coupling. This type of research will bring us closer to understanding invasive propertiesof cells that lead to metastatic cancer and hopefully bring us closer to designing optimal patient-specific treatments.124Chapter 9Concluding remarks9.1 Research summaryThe overarching theme in this research has been the coupling between different levels of organiza-tion that facilitate cell migration. I modelled and analyzed migration patterns in systems rangingfrom developmental biology to cancer, showing how individual cell behaviour can be inferred fromtissue-level observations and vice versa.• By analyzing a simple model, where physical cell properties have been neglected, I demon-strated that drift (such as chemotaxis), could disrupt the formation of cellular aggregates. Theanalysis also provides a link between micro and macro-scale parameters.• During zebrafish development, a cluster of cells (PLLP) needs to segregate and migrate overa long distance. In a model of the proposed interactions between two chemical networks, Wntand FGF, I showed that given an initial asymmetry (i.e. in the steady-state concentration ofWnt receptors) front and back phenotypes could emerge. Coupling this biochemical modelin the HyDiCell3D framework (with deformable cells that migrate based on adhesion andchemotaxis forces) suggested that for long range directional migration, a non-uniform degra-dation of an extracellular chemokine (CXCL12a) was essential. Furthermore, a balance ofcell-cell interactions and chemotaxis forces was required for coordinated cohesive migrationof the PLLP. Model predictions were validated in the lab of Ajay Chitnis (National Instituteof Health).• In modelling collective migration of cancer cells, I found that cells with distinct adhesionproperties can take on leader cell behaviour. The modelling framework guided the design ofexperiments being conducted in the Roskelley lab (UBC). Furthermore, a phenomenologicalmodel of the location of a cell-surface molecule, podocalyxin, coupled to the polarity of a cellcould explain the various cluster morphologies from experiments with distinct types of ECM.125The results suggested that proliferation (specifically, the axis of cell division) influences themorphologies. Lastly, to compare modelling results and biological data, I developed quanti-tative methods for image analysis and automated cell tracking.9.2 ApplicationsFeatures of cell migration can be compared across diverse systems. They are involved in develop-mental biology, normal cell functions (such as wound healing and immune response) and diseases(including cancer). Here, I have studied cell migration in the context of zebrafish development,cell-immune system interactions, a wound healing assay and collective migration of cancer cells.Obviously there are some key differences between these systems, but I found that in the context ofcell migration there are many similarities, including:X Cells with distinct properties develop within a collection of cells, for example the emergenceof leader-cells during wound healing and the segmentation of a cell cluster in zebrafish devel-opment.X The same signalling pathways (CXCL12, Wnt and FGF) that lead to cluster polarization inzebrafish development have been shown to play a role in cancer progression, migration andmetastasis.X The mechanical properties, such as cell adhesion and polarity, are comparable.9.3 Future workThe 3D computational framework, HyDiCell3D, although a good tool to address questions regard-ing cell-cell interactions in a multi-cellular system, is quite complex and not publicly available.Therefore, it may be worth comparing it to a simplified version of cell-cell interactions where anumber of model variants could be explored more easily. Furthermore, this framework would bene-fit from a systematic comparison to other, publicly available, cell-based computational frameworksInteractions between cells and the micro-environment play a key role in facilitating cell migra-tion. In chapter 8, I indicated a number of experiments that could assist our future understanding ofECM mediated cancer cell invasion. Furthermore, the alignment and density of structural fibers inthe ECM have not been captured in my models but it would be an interesting avenue to explore infuture work ( i.e. suggested improvement in Fig. 8.19).Finally, in this thesis (chapter 5) and in my previous research I have studied interactions betweencancer cells and the immune system. I would like to develop a dynamical model that incorporatesthese interactions and uses genomics data to guide the biochemical and biophysical properties ofthe cells. The model would be designed to predict treatment efficacy and plan treatment protocols126for cancer immunotherapies. These treatments have proven successful in a subset of cancers but theresults vary widely and I believe the field would greatly benefit from a mechanistic computationalmodel that can integrate complex interactions.127Bibliography[1] R. Allena and P. Maini. Reaction–diffusion finite element model of lateral line primordiummigration to explore cell leadership. Bulletin of Mathematical Biology, 76(12):3028–3050,2014. → pages 55, 79, 80[2] A. Aman and T. Piotrowski. Wnt/β -Catenin and FGF Signaling Control Collective CellMigration by Restricting Chemokine Receptor Expression. Dev. Cell, 15(5):749–761, 2008.→ pages 7, 8, 54, 56, 59, 71, 76, 77[3] A. Aman, M. Nguyen, and T. Piotrowski. Wnt/β -catenin dependent cell proliferationunderlies segmented lateral line morphogenesis. Developmental Biology, 349(2):470–482,2011. → pages 71, 74[4] R. E. Baker and M. J. Simpson. Models of collective cell motion for cell populations withdifferent aspect ratio: diffusion, proliferation and travelling waves. Physica A: StatisticalMechanics and its Applications, 391(14):3729–3750, 2012. → pages 18[5] A. Beenken and M. Mohammadi. The FGF family: biology, pathophysiology and therapy.Nat. Rev. Drug Discov., 8(3):235–253, 2009. → pages 61, 144[6] D. M. Bryant, J. Roignot, A. Datta, A. W. Overeem, M. Kim, W. Yu, X. Peng, D. J.Eastburn, A. J. Ewald, Z. Werb, et al. A molecular switch for the orientation of epithelialcell polarization. Developmental cell, 31(2):171–187, 2014. → pages 9, 98, 99, 100, 102,117, 119, 120[7] S. W. Byers, C. L. Sommers, B. Hoxter, A. M. Mercurio, and A. Tozeren. Role ofe-cadherin in the response of tumor cell aggregates to lymphatic, venous and arterial flow:measurement of cell-cell adhesion strength. Journal of cell science, 108(5):2053–2064,1995. → pages 107[8] A. E. Carpenter, T. R. Jones, M. R. Lamprecht, C. Clarke, I. H. Kang, O. Friman, D. A.Guertin, J. H. Chang, R. A. Lindquist, J. Moffat, et al. Cellprofiler: image analysis softwarefor identifying and quantifying cell phenotypes. Genome biology, 7(10):R100, 2006. →pages 93[9] P.-C. Chavy-Waddy and T. Kolokolnikov. A local pde model of aggregation formation inbacterial colonies. Nonlinearity, 29(10):3174, 2016. → pages 19, 20, 21, 22, 24, 25, 29, 32,33, 34, 38, 39, 43, 45, 53, 140128[10] K. J. Cheung, E. Gabrielson, Z. Werb, and A. J. Ewald. Collective invasion in breast cancerrequires a conserved basal epithelial program. Cell, 155(7):1639–1651, 2013. → pages 86[11] K. J. Cheung, V. Padmanaban, V. Silvestri, K. Schipper, J. D. Cohen, A. N. Fairchild, M. A.Gorin, J. E. Verdone, K. J. Pienta, J. S. Bader, et al. Polyclonal breast cancer metastasesarise from collective dissemination of keratin 14-expressing tumor cell clusters.Proceedings of the National Academy of Sciences, 113(7):E854–E863, 2016. → pages 9[12] A. Chitnis and D. Dalle Nogare. A computational model reveals the remarkable patterningpotential of the wnt–fgf gene regulatory network in the posterior lateral line primordium.Developmental Biology, 356(1):111, 2011. → pages 8, 55, 72, 79, 80[13] A. B. Chitnis, D. Dalle Nogare, and M. Matsuda. Building the posterior lateral line systemin zebrafish. Developmental Neurobiology, 72(3):234–255, 2012. → pages 56[14] D. Choquet, D. P. Felsenfeld, and M. P. Sheetz. Extracellular matrix rigidity causesstrengthening of integrin–cytoskeleton linkages. Cell, 88(1):39–48, 1997. → pages 72[15] J. Condeelis and J. W. Pollard. Macrophages: obligate partners for tumor cell migration,invasion, and metastasis. Cell, 124(2):263–6, Jan 2006. doi:10.1016/j.cell.2006.01.007. →pages 10[16] D. Dalle Nogare, K. Somers, S. Rao, M. Matsuda, M. Reichman-Fried, E. Raz, and A. B.Chitnis. Leading and trailing cells cooperate in collective migration of the zebrafishposterior lateral line primordium. Development, 141(16):3188–3196, 2014. → pages 8, 55,74, 75, 76, 79, 80, 143[17] C. Dambly-Chaudie`re, N. Cubedo, and A. Ghysen. Control of cell migration in thedevelopment of the posterior lateral line: antagonistic interactions between the chemokinereceptors cxcr4 and cxcr7/rdc1. BMC Developmental Biology, 7(1):23, 2007. → pages 8,54, 72[18] U. Dammer, O. Popescu, P. Wagner, D. Anselmetti, H.-J. Gu¨ntherodt, and G. N. Misevic.Binding strength between cell adhesion proteoglycans measured by atomic forcemicroscopy. Science, 267(5201), 1995. → pages 107[19] J. Delile, M. Herrmann, N. Peyrie´ras, and R. Doursat. A cell-based computational model ofearly embryogenesis coupling mechanical behaviour and gene regulation. Naturecommunications, 8, 2017. → pages 6[20] E. Di Costanzo, R. Natalini, and L. Preziosi. A hybrid mathematical model forself-organizing cell migration in the zebrafish lateral line. Journal of Mathematical Biology,71(1):171–214, 2015. → pages 55, 73, 79, 80, 81[21] D. Drasdo. Center-based Single-cell Models: An Approach to Multi-cellular OrganizationBased on a Conceptual Analogy to Colloidal Particles, pages 171–196. Birkha¨user Basel,Basel, 2007. ISBN 978-3-7643-8123-3. doi:10.1007/978-3-7643-8123-3 8. URLhttps://doi.org/10.1007/978-3-7643-8123-3 8. → pages 6129[22] O. Du Roure, A. Saez, A. Buguin, R. H. Austin, P. Chavrier, P. Siberzan, and B. Ladoux.Force mapping in epithelial cell migration. Proceedings of the National Academy ofSciences of the United States of America, 102(7):2390–2395, 2005. → pages 72[23] L. Edelstein-Keshet and G. B. Ermentrout. Models for spatial polymerization dynamics ofrod-like polymers. Journal of mathematical biology, 40(1):64–96, 2000. → pages 80[24] L. Edelstein-Keshet, W. R. Holmes, M. Zajac, and M. Dutot. From simple to detailedmodels for cell polarization. Philosophical Transactions of the Royal Society of London B:Biological Sciences, 368(1629):20130003, 2013. → pages 101[25] G. B. Ermentrout and L. Edelstein-Keshet. Cellular automata approaches to biologicalmodeling. Journal of theoretical Biology, 160(1):97–133, 1993. → pages 5[26] A. J. Ewald, A. Brenot, M. Duong, B. S. Chan, and Z. Werb. Collective epithelial migrationand cell rearrangements drive mammary branching morphogenesis. Developmental cell, 14(4):570–581, 2008. → pages 9[27] R. J. Filion and A. S. Popel. Intracoronary administration of fgf-2: a computational modelof myocardial deposition and retention. American Journal of Physiology-Heart andCirculatory Physiology, 288(1):H263–H279, 2005. → pages 143[28] A. G. Fletcher, M. Osterfield, R. E. Baker, and S. Y. Shvartsman. Vertex models of epithelialmorphogenesis. Biophysical journal, 106(11):2291–2304, 2014. → pages 6[29] M. J. Footer, J. W. Kerssemakers, J. A. Theriot, and M. Dogterom. Direct measurement offorce generation by actin filament polymerization using an optical trap. Proceedings of theNational Academy of Sciences, 104(7):2181–2186, 2007. → pages 72[30] C. Foroni, M. Broggini, D. Generali, and G. Damia. Epithelial–mesenchymal transition andbreast cancer: Role, molecular mechanisms and clinical impact. Cancer treatment reviews,38(6):689–697, 2012. → pages 8[31] A. Galante and D. Levy. Modeling selective local interactions with memory. Physica D:Nonlinear Phenomena, 260:176–190, 2013. → pages 19[32] T. S. Gardner, C. R. Cantor, and J. J. Collins. Construction of a genetic toggle switch inescherichia coli. Nature, 403(6767):339–342, 2000. → pages 59[33] P. Gerlee and A. R. Anderson. An evolutionary hybrid cellular automaton model of solidtumour growth. Journal of theoretical biology, 246(4):583–603, 2007. → pages 5[34] N. Gompel, N. Cubedo, C. Thisse, B. Thisse, C. Dambly-Chaudie`re, and A. Ghysen. Patternformation in the lateral line of zebrafish. Mechanisms of Development, 105(1–2):69 – 77,2001. → pages 7, 54[35] S. Goswami, E. Sahai, J. B. Wyckoff, M. Cammer, D. Cox, F. J. Pixley, E. R. Stanley, J. E.Segall, and J. S. Condeelis. Macrophages promote the invasion of breast carcinoma cells viaa colony-stimulating factor-1/epidermal growth factor paracrine loop. Cancer Res, 65(12):5278–83, Jun 2005. doi:10.1158/0008-5472.CAN-04-1853. → pages 18, 48, 53130[36] M. L. Graves. Modeling mammary epithelial cell polarization and the role of podocalyxinin breast tumor progression. PhD thesis, University of British Columbia, 2008. URLhttps://open.library.ubc.ca/cIRcle/collections/24/items/1.0066552. → pages 102[37] M. L. Graves, J. A. Cipollone, P. Austin, E. M. Bell, J. S. Nielsen, C. B. Gilks, K. M.McNagny, and C. D. Roskelley. The cell surface mucin podocalyxin regulates collectivebreast tumor budding. Breast Cancer Research, 18(1):11, 2016. → pages 9, 10, 113, 118,119, 122[38] P. Haas and D. Gilmour. Chemokine signaling mediates self-organizing tissue migration inthe zebrafish lateral line. Developmental Cell, 10(5):673–680, 2006. → pages 8, 54, 72, 143[39] C. Hidalgo-Carcedo, S. Hooper, S. I. Chaudhry, P. Williamson, K. Harrington, B. Leitinger,and E. Sahai. Collective cell migration requires suppression of actomyosin at cell–cellcontacts mediated by ddr1 and the cell polarity regulators par3 and par6. Nature cellbiology, 13(1), 2011. → pages 7[40] T. Hillen and K. Painter. Global existence for a parabolic chemotaxis model with preventionof overcrowding. Advances in Applied Mathematics, 26(4):280–301, 2001. → pages 38[41] T. Hillen and K. J. Painter. A user’s guide to pde models for chemotaxis. Journal ofmathematical biology, 58(1-2):183–217, 2009. → pages 18, 19, 38[42] W. R. Holmes and L. Edelstein-Keshet. Analysis of a minimal rho-gtpase circuit regulatingcell shape. Physical biology, 13(4):046001, 2016. → pages 101, 102, 103[43] O. Ilina and P. Friedl. Mechanisms of collective cell migration at a glance. Journal of cellscience, 122(18):3203–3208, 2009. → pages 9[44] B. N. Kholodenko. Cell-signalling dynamics in time and space. Nature reviews MolecularCell Biology, 7(3):165–176, 2006. → pages 56[45] B. Kirkpatrick, L. Nguyen, G. Kondrikova, S. Herberg, and W. D. Hill. Stability of humanstromal-derived factor-1α (CXCL12α) after blood sampling. Ann. Clin. Lab. Sci., 40(3):257–260, 2010. → pages 61, 144[46] H. Knutsdottir, J. S. Condeelis, and E. Palsson. 3-d individual cell based computationalmodeling of tumor cell–macrophage paracrine signaling mediated by egf and csf-1gradients. Integrative Biology, 2016. → pages 137[47] H. Knutsdottir, C. Zmurchok, D. Bhaskar, E. Palsson, D. Dalle Nogare, A. Chitnis, andL. Edelstein-Keshet. Polarization and migration in the zebrafish posterior lateral linesystem. PLoS computational biology, 13(4):e1005451, 2017. → pages 57, 66[48] Y. Kogan, K. E. Halevi-Tobias, G. Hochman, A. K. Baczmanska, L. Leyns, and Z. Agur. Anew validated mathematical model of the wnt signalling pathway predicts effectivecombinational therapy by sfrp and dkk. Biochem. J., 444(1):115–125, 2012. → pages 61,144131[49] V. Lecaudey, G. Cakan-Akdogan, W. H. Norton, and D. Gilmour. Dynamic FGF signalingcouples morphogenesis and migration in the zebrafish lateral line primordium.Development, 135(16):2695–2705, 2008. → pages 7, 54[50] E. Leung, A. Xue, Y. Wang, P. Rougerie, V. Sharma, R. Eddy, D. Cox, and J. Condeelis.Blood vessel endothelium-directed tumor cell streaming in breast tumors requires thehgf/c-met signaling pathway. Oncogene, 36(19):2680, 2017. → pages 10, 18[51] M. M. Lotz, C. A. Burdsal, H. P. Erickson, and D. R. McClay. Cell adhesion to fibronectinand tenascin: quantitative measurements of initial binding and subsequent strengtheningresponse. The Journal of Cell Biology, 109(4):1795–1805, 1989. → pages 107[52] K. E. Luker, J. M. Steele, L. A. Mihalko, P. Ray, and G. D. Luker. Constitutive andchemokine-dependent internalization and recycling of cxcr7 in breast cancer cells todegrade chemokine ligands. Oncogene, 29(32):4599–4610, 2010. → pages 72[53] A. F. Mare´e, V. A. Grieneisen, and P. Hogeweg. The cellular potts model and biophysicalproperties of cells, tissues and morphogenesis. In Single-cell-based models in biology andmedicine, pages 107–136. Springer, 2007. → pages 5[54] A. F. Mare´e, V. A. Grieneisen, and L. Edelstein-Keshet. How cells integrate complexstimuli: the effect of feedback from phosphoinositides and cell shape on cell polarizationand motility. PLoS computational biology, 8(3):e1002402, 2012. → pages 101[55] S. Mark, R. Shlomovitz, N. S. Gov, M. Poujade, E. Grasland-Mongrain, and P. Silberzan.Physical model of the dynamic instability in an expanding cell culture. Biophysical journal,98(3):361–370, 2010. → pages 96[56] M. Matsuda and A. B. Chitnis. Atoh1a expression must be restricted by notch signaling foreffective morphogenesis of the posterior lateral line primordium in zebrafish. Development,137(20):3477–3487, 2010. → pages 143[57] M. Matsuda, D. D. Nogare, K. Somers, K. Martin, C. Wang, and A. B. Chitnis. Lef1regulates dusp6 to influence neuromast formation and spacing in the zebrafish posteriorlateral line primordium. Development, 140(11):2387–97, 2013. → pages 74[58] R. Mayor and S. Etienne-Manneville. The front and rear of collective cell migration. NatureReviews Molecular Cell Biology, 17(2):97–109, 2016. → pages 83[59] J. McGrath, Y. Tardy, C. Dewey, J. Meister, and J. Hartwig. Simultaneous measurements ofactin filament turnover, filament fraction, and monomer diffusion in endothelial cells.Biophysical Journal, 75(4):2070–2078, 1998. → pages 144[60] G. R. Mirams, C. J. Arthurs, M. O. Bernabeu, R. Bordas, J. Cooper, A. Corrias, Y. Davit,S.-J. Dunn, A. G. Fletcher, D. G. Harvey, et al. Chaste: an open source c++ library forcomputational physiology and biology. PLoS computational biology, 9(3):e1002970, 2013.→ pages 6132[61] P. J. Murray, C. M. Edwards, M. J. Tindall, and P. K. Maini. From a discrete to a continuummodel of cell dynamics in one dimension. Physical Review E, 80(3):031912, 2009. →pages 18, 19[62] K. Nabeshima, T. Inoue, Y. Shimao, and T. Sameshima. Matrix metalloproteinases in tumorinvasion: role for cell migration. Pathology international, 52(4):255–264, 2002. → pages120[63] A. Nechiporuk and D. W. Raible. Fgf-dependent mechanosensory organ patterning inzebrafish. Science, 320(5884):1774–1777, 2008. → pages 7, 54[64] J. S. Nielsen and K. M. McNagny. The role of podocalyxin in health and disease. Journal ofthe American Society of Nephrology, 20(8):1669–1676, 2009. → pages 9[65] M. A. Nugent and E. R. Edelman. Kinetics of basic fibroblast growth factor binding to itsreceptor and heparan sulfate proteoglycan: a mechanism for cooperativity. Biochemistry, 31(37):8876–8883, 1992. → pages 61, 143, 144[66] H. G. Othmer, S. R. Dunbar, and W. Alt. Models of dispersal in biological systems. Journalof mathematical biology, 26(3):263–298, 1988. → pages 19[67] A. W. Overeem, D. M. Bryant, and S. C. van IJzendoorn. Mechanisms of apical–basal axisorientation and epithelial lumen positioning. Trends in cell biology, 25(8):476–485, 2015.→ pages 122[68] K. J. Painter. Continuous models for cell migration in tissues and applications to cell sortingvia differential chemotaxis. Bulletin of mathematical biology, 71(5):1117–1147, 2009. →pages 38[69] K. J. Painter, N. J. Armstrong, and J. A. Sherratt. The impact of adhesion on cellularinvasion processes in cancer and development. Journal of theoretical biology, 264(3):1057–1067, 2010. → pages 38[70] E. Palsson. A 3-d model used to explore how cell adhesion and stiffness affect cell sortingand movement in multicellular systems. Journal of Theoretical Biology, 254(1):1–13, 2008.→ pages 6, 12, 16, 73, 138[71] E. Palsson and H. G. Othmer. A model for individual and collective cell movement indictyostelium discoideum. Proceedings of the National Academy of Sciences, 97(19):10448–10453, 2000. ISSN 0027-8424. doi:10.1073/pnas.97.19.10448. URLhttp://www.pnas.org/content/97/19/10448. → pages 6[72] P. Panorchan, M. S. Thompson, K. J. Davis, Y. Tseng, K. Konstantopoulos, and D. Wirtz.Single-molecule analysis of cadherin-mediated cell-cell adhesion. Journal of cell science,119(1):66–74, 2006. → pages 91[73] J. Panovska-Griffiths, K. M. Page, and J. Briscoe. A gene regulatory motif that generatesoscillatory or multiway switch outputs. Journal of The Royal Society Interface, 10(79):20120826, 2013. → pages 59, 64133[74] M. J. Paszek, N. Zahir, K. R. Johnson, J. N. Lakins, G. I. Rozenberg, A. Gefen, C. A.Reinhart-King, S. S. Margulies, M. Dembo, D. Boettiger, et al. Tensional homeostasis andthe malignant phenotype. Cancer cell, 8(3):241–254, 2005. → pages 98, 122[75] A. Patsialou, J. Wyckoff, Y. Wang, S. Goswami, E. R. Stanley, and J. S. Condeelis. Invasionof human breast cancer cells in vivo requires both paracrine and autocrine loops involvingthe colony-stimulating factor-1 receptor. Cancer Res, 69(24):9498–506, Dec 2009.doi:10.1158/0008-5472.CAN-09-1868. → pages 10, 48, 53[76] L. Petitjean, M. Reffay, E. Grasland-Mongrain, M. Poujade, B. Ladoux, A. Buguin, andP. Silberzan. Velocity fields in a collectively migrating epithelium. Biophysical journal, 98(9):1790–1800, 2010. → pages 9[77] I. Ramis-Conde, D. Drasdo, A. R. Anderson, and M. A. Chaplain. Modeling the influenceof the e-cadherin-β -catenin pathway in cancer cell invasion: a multiscale approach.Biophysical journal, 95(1):155–165, 2008. → pages 6[78] M. Reffay, M.-C. Parrini, O. Cochet-Escartin, B. Ladoux, A. Buguin, S. Coscoy, F. Amblard,J. Camonis, and P. Silberzan. Interplay of rhoa and mechanical forces in collective cellmigration driven by leader cells. Nature cell biology, 16(3):217–223, 2014. → pages 100[79] L. G. Reyna and M. J. Ward. Metastable internal layer dynamics for the viscouscahn-hilliard equation. Methods and Applications of Analysis, 2:285–306, 1995. → pages29[80] R. Riahi, Y. Yang, D. D. Zhang, and P. K. Wong. Advances in wound-healing assays forprobing collective cell migration. Journal of laboratory automation, 17(1):59–65, 2012. →pages 9[81] P. G. Saffman and F. S. G. TAYLOR. The penetration of a fluid into a porous medium orhele-shaw cell containing a more viscous liquid. In Dynamics of Curved Fronts, pages155–174. Elsevier, 1988. → pages 96[82] Y. Saka and J. C. Smith. A mechanism for the sharp transition of morphogen gradientinterpretation in xenopus. BMC Developmental Biology, 7(1):47, 2007. → pages 56, 59, 65[83] M. Scianna and L. Preziosi. A hybrid model describing different morphologies of tumorinvasion fronts. Mathematical modelling of natural phenomena, 7(1):78–104, 2012. →pages 86[84] M. Scianna, L. Preziosi, and K. Wolf. A cellular potts model simulating cell migration onand in matrix environments. Mathematical Biosciences and Engineering, 10(1):235–261,2013. → pages 5[85] E. R. Shamir and A. J. Ewald. Three-dimensional organotypic culture: experimental modelsof mammalian biology and disease. Nature reviews Molecular cell biology, 15(10):647–664, 2014. → pages 1134[86] A. Somasiri, J. S. Nielsen, N. Makretsov, M. L. McCoy, L. Prentice, C. B. Gilks, S. K. Chia,K. A. Gelmon, D. B. Kershaw, D. G. Huntsman, et al. Overexpression of the anti-adhesinpodocalyxin is an independent predictor of breast cancer progression. Cancer research, 64(15):5068–5073, 2004. → pages 10, 98[87] M. S. Steinberg and M. Takeichi. Experimental specification of cell sorting, tissuespreading, and specific spatial patterning by quantitative differences in cadherin expression.Proceedings of the National Academy of Sciences, 91(1):206–209, 1994. → pages 88[88] S. J. Streichan, G. Valentin, D. Gilmour, and L. Hufnagel. Collective cell migration guidedby dynamically maintained gradients. Physical Biology, 8(4):045004, 2011. → pages 8, 55,79[89] M. H. Swat, G. L. Thomas, J. M. Belmonte, A. Shirinifard, D. Hmeljak, and J. A. Glazier.Multi-scale modeling of tissues using compucell3d. Methods in cell biology, 110:325, 2012.→ pages 6[90] T. Takeda, W. Y. Go, R. A. Orlando, and M. G. Farquhar. Expression of podocalyxininhibits cell–cell adhesion and modifies junctional properties in madin-darby canine kidneycells. Molecular biology of the cell, 11(9):3219–3232, 2000. → pages 107[91] G. I. Taylor et al. Dispersion of soluble matter in solvent flowing slowly through a tube.Proc. R. Soc. Lond. A, 219(1137):186–203, 1953. → pages 37[92] E. Theveneau, B. Steventon, E. Scarpa, S. Garcia, X. Trepat, A. Streit, and R. Mayor.Chase-and-run between adjacent cell populations promotes directional collective migration.Nature cell biology, 15(7):763, 2013. → pages 51[93] K. Twaroski, S. K. Mallanna, R. Jing, F. DiFurio, A. Urick, and S. A. Duncan. Fgf2mediates hepatic progenitor cell formation during human pluripotent stem celldifferentiation by inducing the wnt antagonist nkd1. Genes & development, 29(23):2463–2474, 2015. → pages 59, 78[94] L. E. Valdivia, R. M. Young, T. A. Hawkins, H. L. Stickney, F. Cavodeassi, Q. Schwarz,L. M. Pullin, R. Villegas, E. Moro, F. Argenton, et al. Lef1-dependent wnt/β -cateninsignalling drives the proliferative engine that maintains tissue homeostasis during lateralline development. Development, 138(18):3931–3941, 2011. → pages 74[95] G. Valentin, P. Haas, and D. Gilmour. The chemokine sdf1a coordinates tissue migrationthrough the spatially restricted activation of cxcr7 and cxcr4b. Current Biology, 17(12):1026–1031, 2007. → pages 8, 54, 72[96] M. Venero Galanternik, Marina, K. L. Kramer, and T. Piotrowski. Heparan sulfateproteoglycans regulate fgf signaling and cell polarity during collective cell migration. CellRep., 10(3):414–428, 2015. → pages 145[97] G. Venkiteswaran, S. W. Lewellis, J. Wang, E. Reynolds, C. Nicholson, and H. Knaut.Generation and dynamics of an endogenous, self-generated signaling gradient across amigrating tissue. Cell, 155(3):674–687, 2013. → pages 72135[98] U. Wilensky and M. Resnick. Thinking in levels: A dynamic systems approach to makingsense of the world. Journal of Science Education and technology, 8(1):3–19, 1999. →pages 6[99] K. Wolf, Y. I. Wu, Y. Liu, J. Geiger, E. Tam, C. Overall, M. S. Stack, and P. Friedl.Multi-step pericellular proteolysis controls the transition from individual to collectivecancer cell invasion. Nature cell biology, 9(8):893–904, 2007. → pages 9[100] J. Wyckoff, W. Wang, E. Y. Lin, Y. Wang, F. Pixley, E. R. Stanley, T. Graf, J. W. Pollard,J. Segall, and J. Condeelis. A paracrine loop between tumor cells and macrophages isrequired for tumor cell migration in mammary tumors. Cancer Res, 64(19):7022–9, Oct2004. doi:10.1158/0008-5472.CAN-04-1449. → pages 18, 48, 53[101] J. B. Wyckoff, Y. Wang, E. Y. Lin, J.-f. Li, S. Goswami, E. R. Stanley, J. E. Segall, J. W.Pollard, and J. Condeelis. Direct visualization of macrophage-assisted tumor cellintravasation in mammary tumors. Cancer research, 67(6):2649–2656, 2007. → pages 10[102] N. Yamaguchi, T. Mizutani, K. Kawabata, and H. Haga. Leader cells regulate collective cellmigration via rac activation in the downstream signaling of integrin β1 and pi3k. Scientificreports, 5, 2015. → pages 9, 83[103] W. Yu, A. Datta, P. Leroy, L. E. O’Brien, G. Mak, T.-S. Jou, K. S. Matlin, K. E. Mostov, andM. M. Zegers. β1-integrin orients epithelial polarity via rac1 and laminin. Molecularbiology of the cell, 16(2):433–445, 2005. → pages 98, 102, 114[104] W. Yu, A. M. Shewan, P. Brakeman, D. J. Eastburn, A. Datta, D. M. Bryant, Q.-W. Fan,W. A. Weiss, M. M. Zegers, and K. E. Mostov. Involvement of rhoa, rock i and myosin ii ininverted orientation of epithelial polarity. EMBO reports, 9(9):923–929, 2008. → pages 118[105] S. Zhou and M. Y. Wang. Multimaterial structural topology optimization with a generalizedcahn–hilliard model of multiphase transition. Structural and MultidisciplinaryOptimization, 33(2):89, 2007. → pages 9136Appendix ASupporting InformationA.1 The computational frameworkA.1.1 Chemical concentrationThe surfaces of the simulated cells act as sources and sinks for ligands (signalling chemicals). Totrack the ligand concentration, the following reaction-diffusion PDE is solved for each ligand, L, ona 3D cartesian grid (with each grid cube representing 103 µm3):∂L∂ t= DL∇2L+ pLΩC−δLL−κLLΩC. (A.1)DL is the diffusion coefficient for ligand L, pL is its rate of secretion by cell C, δL is the uniformglobal decay, κL is the local degradation of ligand L by cell C and ΩC represents the proportion ofcell surface area in each lattice cube (0≤ΩC ≤ 1).Cells are not point sources but instead can overlap several lattice cubes in the simulations. Thisis accounted for in the model and the details can be found in [46]. However, in the setup of themodels considered in this thesis, this may not greatly affect the results.A.1.2 ForcesThe adhesion force between two cells acts along the same axis, with a magnitude that goes to zerofor cells that overlap or that are far apart. This force pulls cells towards their equilibrium spacing,137but acts locally. The net force at cell interfaces, Finter f ace, is calculated as:Finter f ace = −Fadhesionχ[(x+ x0)e−λ (x+x0)2)− v0e−λx2]· ~ri j‖~ri j‖ if x > 0,Fexclusionχ(−x) 95 · ~ri j‖~ri j‖ if x≤ 0.χ =rcell2(1di+1d j), x =d−dminrcell, (A.2)where~ri j is a vector between the centers of cells i and j, rcell is the cell radius, di and d j are thedistances from the centre to the cell surface along the vector~ri j for ellipsoid i and j, respectively,d is the distance between the surfaces of two neighboring ellipsoids along the vector~ri j (i.e. d =‖~ri j‖−di−d j) and dmin determines the point at which the cells are in the equilibrium state where theadhesion and exclusion terms balance (see Fig. 11 in [70]). In all simulations dmin =−0.1 ·rcell . Thisallows the ellipsoids to overlap and corrects for the fact that cells are not actually rigid ellipsoids.(It also allows neighboring cells to have a finite contact area, which would not be the case forimpermeable rigid ellipoids.) The constants x0 and v0 are chosen to ensure that when cells aretouching, i.e. when x = 0, both Fadhesive = 0 and F ′adhesive(x) = 0. The interface force, Finter f ace, isplotted in Fig. 3.2.A) B)Figure A.1: A) illustration of the viscolesticity of the ellipsoid obtained using a Maxwell el-ement (spring, with spring constant κ1, in series with a dashpot whose frictional coef-ficient is µ1) in parallel with a spring (with spring constant κ2). B) illustration of thedistances between two cells (i and j) used to calculate the forces in Eqn. (A.2)A.1.3 Cell shapeThe cell is modelled as a viscoelastic ellipsoid with semi-minor axes a,b and c. Each axis containsa Maxwell element (spring, with spring constant κ1, in series with a dashpot whose frictional co-138efficient is µ1) in parallel with a spring (with spring constant κ2). (See [24] for typical mechanicalparameters.) The forces acting on a cell are resolved onto the three ellipsoidal axes and the shapeof the cell is obtained from these mechanical elements under volume conservation.dridt=κ1(Fi−Fmod)µ1(κ1+κ2)+dFi/dt(κ1+κ2)− ri κ1κ2µ1(κ1+κ2) , (A.3)rarbrc = (ra+∆ra)(rb+∆rb)(rc+∆rc) =Vellipse/(43pi), (A.4)where ri is the length of either the a, b or c axis in units of 10 µm and Fmod is calculated fromthe volume constraint in each time step by solving Eqn. (A.3) to find ∆ra, ∆rb and ∆rc under theconstraint of Eqn. (A.4). Allowing cells to deform is significant, as it is easier for deformable cellsto move past neighbors in 3D systems.A.1.4 Cell divisionCell growth is modelled with growth of each of the three primary axes of the ellipsoid. Furthermore,if a cell is in a tightly packed environment it experiences pressure from its neighbours. In themodel, this pressure decreases the rate at which the cell grows. More specifically, the growth rateis proportional to 1010+P2 , where P is the pressure force. A cell will divide when the volume hasdoubled.A.2 Additional calculations for analysis in chapter 5Rescaled equation for two interacting cell typesI rescale space in Eqn. (5.17) as x=α1xˆ. I then drop the hat notation and get:ut = uxx−B1(wxxuw− (wx)2 uw2 +wxuxw),wt = Awxx−B2(uxxwu− (ux)2 wu2 +uxwxu),(A.5)whereA =α2α1=2a2+ c21− c22d2a1+ c12− c11d ,B1 =β1α1=c12(d+1)2a1+ c12− c11d , andB2 =β2α1=c21(d+1)2a1+ c12− c11d .Inhomogeneous steady stateFollowing is the O(h4) terms for the continuum limit of the model in chapter 5:139g(iv)(0) =−uwxxxxw(d+1+d(d+1)22)+wxxxwxuw2(4(d+1)+d(d+1)22+(d+1)2(2d+1)2)− wxx(wx)2uw3(12(d+1)+(d+1)2(2d+1))+(wxx)2uw2(3(d+1)+(d+1)2(2d+1)2)+(wx)4uw46(d+1)− uxwxxxw(3(d+1)+d(d+1)22)+uxwxwxxw2(9(d+1)+(d+1)2(2d+1)2)− ux(wx)3w36(d+1)+uxx(wx)2w23(d+1)− uxxwxxw3(d+1)− uxxxwxw(d+1)+uxxxx2(A.6)A quick check shows that when u = w , the system from [9] is retrieved.The continuum limit for U in Eqn. (5.1) up to O(h4) is:ut = a1(h2uxx+h412uxxxx)+ c11(−h2d2uxx+h412([−12− d2(d2+2d+3)]uxxxx+(d+1)2(2d+1)2(uxuxxu)x))+ c12(h2(uxx2− u2wwxx(d+1)+(wx)2u(d+1)2w2−wx ux(d+1)2w)+h412(g(iv)(0))).(A.7)140A.3 Further details for zebrafish researchA.3.1 Model equations, derivation and scalingThe receptor equations, from Eqs. (6.4) aredFRdt= EF(WB)− rFFR = EF,11+(WB/W0)n − rFFR (A.8a)dWRdt= EW (FB)− rWWR = EW,11+(FB/F0)m − rWWR. (A.8b)I will recast this in terms of FR,WR alone. Based on Michaelis-Menten kinetics, Eqs. (6.3), thebound receptors are given byFB(t) =FR(t)F(t)KF +F= FR(t)(F/KF)1+(F/KF)= FR(t)sF(F/KF), (A.9a)WB(t) =WR(t)W (t)KW +W=WR(t)(W/KW )1+(W/KW )=WR(t)sW (W/KW ), (A.9b)where sF and sW are defined as:sF(x) =x1+ x, sW (x) =x1+ x.Now the ligand dependence is in the expressions sF ,sW which represent the fraction of bound re-ceptors for a given level of ligand concentration (relative to the ligand KD). Note that both functionssatisfy 0 ≤ s j(x) ≤ 1. I scale the FGF (respectively Wnt) ligand concentration by KF (respectivelyKW ).Substituting these expressions into the receptor dynamics equations leads todFRdt= rF[(F11+(WRsW/W0)m)−FR], (A.10a)dWRdt= rW[W11+(FRsF/F0)n−WR]. (A.10b)I scale the FGF receptors by F1 and the Wnt receptors by W1. Then the result isdFRdt= rF[11+(WR/ω)m−FR], (A.11a)dWRdt= rW[11+(FR/φ)n−WR], (A.11b)141whereφ = F0/(F1sF), ω =W0/(W1sW ). (A.11c)Note that as ligand concentration F (respectively W ) decreases, sF (respectively sW ) decreases, andso φ (respectively ω) increases. Now Eqs. (A.11) a,b are in the form of a mutual inhibition moduleand have a number of possible behaviours that depend on parameters φ and ω , as shown in Fig. 6.4.In 1D, the ligand equations are∂F∂ t= DF∂ 2F∂x2−δFF− rate of binding + rate of production (A.12a)∂W∂ t= DW∂ 2W∂x2−δWW − rate of binding + rate of production (A.12b)Ligand can bind to unoccupied receptors. Then, based on mass-action, the rates of binding arerate of binding to FGF receptors = kF,onF(FR−FB) = kF,onFFR(1− FKF +F ) (A.12c)rate of binding to Wnt receptors = kW,onW (WR−WB) = kW,onWWR(1− WKW +W ) (A.12d)where the constants KF ,KW are given byKF =kF,off+ kF,2kF,on, KW =kW,off+ kW,2kW,on. (A.12e)See Fig. 6.2 for these rate parameters. The binding rates can be simplified tokF,onFFR(1− FKF +F ) = kF,onFFRKFKF +F= kF,onFFR11+(F/KF), (A.12f)kW,onWWR(1− WKW +W ) = kW,onWWRKWKW +W= kW,onWWR11+(W/KW ). (A.12g)Similarly, the rates of production arerate of production of FGF ligand = pFWB = pFWRWKW +W(A.13)rate of production of Wnt ligand = pWWB = pWWRWKW +W. (A.14)The reaction-diffusion PDEs for ligand concentrations can be written in the form∂F∂ t= DF∂ 2F∂x2−δFF− kF,onFFR 11+(F/KF) + pFWRWKW +W, (A.15a)∂W∂ t= DW∂ 2W∂x2−δWW − kW,onWWR 11+(W/KW ) + pWWRWKW +W. (A.15b)142I now put these in dimensionless form by the following scaling:F = F∗F¯ , W =W ∗W¯ ,where *s are dimensionless and bars are the scales. I choose F¯ =KF ,W¯ =KW as noted earlier. I alsohave scales for distance (x = x∗L) and for the receptor levels (FR = F∗R F1,WR =W∗RW1) as discussedpreviously. This gives the dimensionless (except for time) model:∂F∂ t=DF∂ 2F∂x2−δFF−κFFR F1+F +ρFWRW1+W, (A.16a)∂W∂ t=DW∂ 2W∂x2−δWW −κWWR W1+W +ρWWRW1+W, (A.16b)whereDi =DiL2, κF = kF,onF1, κW = kW,onW1, ρF =pFW1KF, and ρW =pWW1KW. (A.16c)A.3.2 Dose response experiment in zebrafish embryoThese experiments were conducted by Damian Dalle Nogare at the National Institute of Health.CldnB:lynGFP [38] embryos were treated from 32-36hpf with the indicated concentrations of SU5402while shaking and then fixed in 4% PFA overnight. In situ hybridization for lef1 transcript was per-formed as previously described [56]. Following in situ hybridization, embryos were stained withanti-GFP antibody to visualize the CldnB:lynGFP. The lef1 ratio was computed by dividing thelength of the lef1 domain by the total length of the PLLP.A.3.3 Parameter Estimation and ValuesPrimordium size and units I use units of µm for length, and minutes for time. The PLLP is on theorder of 100-200 µm long [16]. The entire primordium is estimated to contain about 125 cells atthe onset of migration.Ligand-receptor binding kinetics and typical concentrations First, I make the assumption that theendocytosis rates of Wnt and FGF receptors are the same: rW = rF . The receptor endocytosis ratesappear in the scaled receptor ODE equations as a time-scaling factor only, and given the lack ofinformation on endocytosis rates, I set rW = rF = 1. I set n = m = 5, but values n 6= m ≥ 2 wouldbe equally suitable to set up the bistable mutual inhibition system.FGF ligand-receptor dynamics have been well studied. For instance, [27, 65] contain detailedinformation in tables of parameter values and detailed assumption for estimated unknown values. In143[65] I find the on and off rates for FGF binding to its cell surface receptor kF,on = 0.227 nM−1min−1.To calculate Michaelis-Menten constant, KF , for FGF binding, the phosphorylation rate of thereceptor-ligand complex needs to be estimated. Using an estimate of 1-2 seconds for receptor-ligandcomplex phosphorylation, the Michaelis-Menten constant isKF =kF,off+ kF,2kF,on=0.003 min−1+0.02 min−10.227 nM−1 min−1≈ 0.1 nM.The decay rate of FGF ligand in the PLLP is estimated using the mean half life for FGF2 in thehuman body, which is 7.6 hours [5]. This estimate givesδF =ln2456 min= 0.0015 min−1.I can find similar information for Wnt ligand-receptor dynamics. For instance, [48] validate amathematical model of the Wnt signalling pathway. In this research, I find the on and off rates forWnt ligand binding to its cell surface receptor, Frizzled, as follows: kW,on = 7.9×104 M−1 s−1 andkW,off = 4.7×10−4 s−1.I use CXCL12a as a proxy for a the unknown Wnt ligand. In [45], the half-life of humanCXCL12a is 26 minutes. This estimate givesδW =ln226 min= 0.027 min−1.The parameters F1 and W1 are the steady-state concentration of FGF and Wnt receptors in theabsence of competition. For these parameters, I expect the number of FGF receptors to be on theorder of thousands per cell. This corresponds to a concentration of approximately 1 nM, hence, Iestimate F1 = 1 nM. From [48], I find that there are 30 Frizzled receptors per cell and set W1 = 0.03nM. The parameters F0 and W0 are the concentrations at which bound FGF receptors and boundWnt receptors reach half-max inhibition of the other type of receptor, respectively. I set W0 = 0.001nM and F0 = 0.1 nM.The Michaelis-Menten constants, KW and KF , can be used to scale the results from my model tocompare to the equivalent biological parameters. However, finding data to estimate these parametershas proven difficult. Since these parameters are scaled out in the analysis and do not affect the resultsI have set them to be 1.Estimating the diffusion rates To estimate the diffusion rates of FGF and Wnt ligand, I compareliterature values and estimate the diffusion rate from molecular weight. FGF is a protein with amolecular weight of 30.7 kD. (By comparison, actin is 42 kD.) Actin diffuses at a rate of roughly5 µm2/s in the cytosol [59], and 10 times faster in water. The ratio of the diffusion coefficientswould be roughly the same as the ratio of cube root of the molecular weights. Using the estimate144of diffusion of actin in water as 50µm2/s and the ratio of cube roots of the molecular weights(42/31)1/3 ≈ 1.1, I find the FGF diffusion coefficient closer to DF ≈ 55µm2/s = 3300µm2/min.However, this estimate seems too fast for FGF ligand diffusion within the PLLP given the size andscale of the PLLP.The binding to receptors and activity of FGF is known to depend on Heparan Sulfate Proteo-glycans (HSPG). Not only have HSPGs been shown to regulate FGF signalling within the PLLP,but also to be an integral part of the feedback loop that organizes the PLLP. In fact, in the PLLPwith mutant HSPGs, wild type signalling is disrupted [96]. Due to the binding of FGF to HSPGs,I consider an effective diffusion coefficient for FGF ligand. Let kon and koff be the forward andbackward rates of FGF binding to HSPGs, DH be the diffusion coefficient of FGF bound to HSPGs,and DF the diffusion coefficient of free FGF.Deffective =DHkon+DFkoffkon+ koff=1τH + τF(DHτH +DFτF), (A.17)where τH = 1koff and τF =1kon. Then the effective diffusion coefficient can be interpreted as theweighted average diffusion where the weights are the mean fractional residence times on the boundto HSPGs and freely diffusing. In this way, the diffusion coefficient of FGF within the PLLP can beregulated by binding to HSPGs. Consequently, I set DF = 10µm2/s = 600µm2/min. Wnt ligandhas a molecular weight similar to actin (41 kD), and so I again use the ratio of cube roots of themolecular weights to estimate that DW ≈ DF = 600µm2/min.145
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The multi-levelled organization of cell migration :...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The multi-levelled organization of cell migration : from individual cells to tissues Knutsdottir, Hildur 2018
pdf
Page Metadata
Item Metadata
Title | The multi-levelled organization of cell migration : from individual cells to tissues |
Creator |
Knutsdottir, Hildur |
Publisher | University of British Columbia |
Date Issued | 2018 |
Description | Cell migration is a complex interplay of biochemical and biophysical mechanisms. I investigate the link between individual and collective cell behaviour using mathematical and computational modelling. Specifically, I study: (1) cell-cell interactions in a discrete framework with a spatial sensing range, (2) migration of a cluster of cells during zebrafish (Danio rerio) development, and (3) collective migration of cancer cells and their interactions with the extracellular-matrix (ECM). My 1D model (1), is approximated by a continuum equation and investigated using asymptotic approximations, steady-state analysis, and linear stability analysis. Analysis and computations characterize regimes corresponding to cell clustering, and provide a link between micro and macro-scale parameters. Results suggest that drift (i.e. due to chemotaxis), can disrupt the formation of cellular aggregates. In (2), I investigate spontaneous polarization of a cell-cluster (the posterior lateral line primordium, PLLP) in zebrafish development. I use a cell-based computational framework (hybrid discrete cell model, HyDiCell3D) coupled with differential equation model to track the segregation and migration of the PLLP. My model includes mutual inhibition between the diffusible growth factors Wnt and FGF. I find that a non-uniform degradation of an extracellular chemokine (CXCL12a) and chemotaxis is essential for long-range cohesive migration. Results compare favourably with data from the Chitnis lab (NIH). I continue using HyDiCell3D in (3) to elucidate mechanisms that facilitate cancer invasion. I focus on: wound healing in a cell-sheet (2D epithelium), and cell-clusters (3D spheroids) embedded in ECM with internal signalling mediated by podocalyxin, a trans-membrane molecule. Experimental data from the Roskelley lab (UBC) motivates the model derivation. I use the models to investigate the role of cell-cell and cell-ECM adhesion in collective migration as well as the emergence of a distinct phenotype (leader-cells) that guide the migration. ECM induced disruption in the localization of podocalyxin on the cell membrane is captured in the model along with morphological changes of spheroids. The model predicts that cell polarity and cell division axis influence the invasive potential. Lastly, I developed quantitative methods for image analysis and automated tracking of cells in a densely packed environment to compare modelling results and biological data. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2018-04-26 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0366016 |
URI | http://hdl.handle.net/2429/65650 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2018-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2018_september_knutsdottir_hildur.pdf [ 59.33MB ]
- Metadata
- JSON: 24-1.0366016.json
- JSON-LD: 24-1.0366016-ld.json
- RDF/XML (Pretty): 24-1.0366016-rdf.xml
- RDF/JSON: 24-1.0366016-rdf.json
- Turtle: 24-1.0366016-turtle.txt
- N-Triples: 24-1.0366016-rdf-ntriples.txt
- Original Record: 24-1.0366016-source.json
- Full Text
- 24-1.0366016-fulltext.txt
- Citation
- 24-1.0366016.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0366016/manifest