University of Amsterdam
Masters Thesis
Influence of fibroblatic reticular cell
network on T cell movement in lymph
nodes
Author:
Laurens Stuurman
Supervisor:
Shabaz Sultan
A thesis submitted in partial fulfilment of the requirements for the degree of Master of Science in Computational Science
in the
Computational Immunology Group Department of Tumor Immunology Radboud University Medical Center
I, Laurens Stuurman, declare that this thesis, entitled ‘Influence of fibroblatic reticular cell network on T cell movement in lymph nodes’ and the work presented in it are my own. I confirm that:
This work was done wholly or mainly while in candidature for a research degree
at the University of Amsterdam.
Where any part of this thesis has previously been submitted for a degree or any
other qualification at this University or any other institution, this has been clearly stated.
Where I have consulted the published work of others, this is always clearly
at-tributed.
Where I have quoted from the work of others, the source is always given. With
the exception of such quotations, this thesis is entirely my own work.
I have acknowledged all main sources of help.
Where the thesis is based on work done by myself jointly with others, I have made
clear exactly what was done by others and what I have contributed myself.
ii
UNIVERSITY OF AMSTERDAM
Abstract
Faculty of Science Informatics Institute
Master of Science in Computational Science
Influence of fibroblatic reticular cell network on T cell movement in lymph nodes
by Laurens Stuurman
One of the first stages of the immune response consists of T cells recognizing foreign antigen. T cells continuously migrate through the body, circulating in blood and lymph nodes in search of antigen presenting dendritic cells. Hence the search behaviour of T cells within lymph nodes is very important for an effective immune response. Movement of T cell is guided by a network of fibroblastic reticular cells (FRCs), but is also influenced by the densely populated environment inside a lymph node. Current project aims to inquire how the structure of the FRC network influences T cell motility.
First part of this project is aimed at generating structurally different models for the FRC network and comparing them with previous research and real data. We find the random geometric graphs produce structures best resembling real FRC networks, both by means of the graph theoretical measures of small-world properties and by quantifying the negative spaces in the network by means of a gap analysis.
The second part of this project is aimed at modelling T cell movement in lymph nodes. To do so, two versions of the Cellular Potts model (CPM) were used; the preferred di-rection model and the act model. Single cell simulations validated that both models can simulate cells moving with similar speed and persistent as real T cells. In a fully popu-lated environment realistic cell movement was harder to simulate. Finally we performed simulations with structurally different FRC networks. We show that lattice like FRC structure allow for more persistent movement then other structures. Also lattice like FRC structures allow for streams of coherently moving T cells to form. We show that streams of coherently moving T cells mainly form along fibers of the network, confirming the hypothesis that the FRC structure forms a ”road map” for T cell migration within the lymph node.
Declaration of Authorship i
Abstract i
Contents iii
List of Figures v
Abbreviations viii
0.1 The adaptive immune response and lymph nodes . . . 1
1 Quantifying and 3D modeling of the fibroblastic reticular cell network structure 4 1.1 Introduction. . . 4
1.2 Experiments. . . 8
1.2.1 Small-world properties in different graph models . . . 8
1.2.2 From graphs to 3D models . . . 8
1.2.3 Gap Analysis . . . 9
1.3 Results. . . 10
1.3.1 Small world properties of different graph models . . . 10
1.3.2 Conclusion . . . 15
2 Simulating migration behaviour of T-cells in lymph node 16 2.1 Introduction. . . 16
2.2 Methods . . . 21
2.2.1 Cellular Potts Model . . . 21
2.2.2 Preferred direction in CPM . . . 21
2.2.3 Act model in CPM . . . 22
2.3 Experiments. . . 23
2.3.1 Single cell simulations . . . 23
2.3.2 Simulations at different densities . . . 23
2.3.3 Simulating fully populated environment . . . 24
2.3.4 Effect of different FRC structures. . . 24
2.3.5 Measures to quantify migration behaviour . . . 25
2.4 Results. . . 28 iii
Contents iv
2.4.1 Single cell simulations . . . 28
2.4.2 Increasing density of simulated environment . . . 30
2.4.3 Simulating fully populated LN . . . 31
2.4.4 Effect of FRC network topology . . . 32
2.5 Discussion and future work . . . 39
1.1 Small-world measures for all types of simulated networks. Observed small-world properties of FRC as observed by Novkovic and colleagues were σ : 6.7 and ω : −0.27 .Points represent average σ and ω values for generated graphs of certain size and density. In the WS and BA-cl models one extra parameter can be varied, which causes the higher amount of data points in these graphs. Size of the points show summed variance within graphs of specific density and size, ranging between .0054 and 1.345. . . 10
1.2 Effect of networks size and density on σ and ω small-world measures. Sizes of the dots indicate standard deviation of σ and ω measures. σ Small-world measure is clearly influenced by both network size and den-sity. Whereas ω measure captures differences between networks more consistently across different network sizes and densities. . . 11
1.3 Degree distributions for generated graphs with similar size and density as FRC network reported in [1] : N = 176 , hki = 8. Specific parameters of graph models : PER = 0.0457, MBA = MBA−cl = 4,PBA−cl = 0.666,
KW S = 8,PW S = 0.0278, rGM −3D= 0.1234. . . 12
1.4 Distribution of edge lengths of 3D models based on different graph types. 13
1.5 Cross-section of stromal structure in mouse LN and simulated stromal cell structure, filled with circles after gap analysis. The yellow structure in B represents the simulated stromal structure.. . . 13
1.6 Proportion of surface area filled with circles of varying radii in real mouse FRC networks and graph based 3D models. Gray lines represent va-riety of gaps present in real mouse LNs. Left pannel shows raw data, right pannel smoothed data using B-spline representation. Error between smoothed curves of grahp models and mouse data was lowest for the ran-dom geometric model (GM-3D). . . 14
2.1 Example of more random like FRC structure (A) and more lattice like FRC structure (B). . . 17
2.2 Difference between only observing a subset of cells, versus the fully pop-ulated environment of the LN. Thin lines in the left panel show represent the track of the center of mass of the simulated cells. Gray structure is the a model of FRC structure based on random geometric graph. Differ-ent colors of cells have no significant meaning. Dimension of simulated environment is 32 cubic microns. . . 18
2.3 FRC models based on different graph types. . . 25
2.4 Average speed and persistence times as resulted from second grid search for single cell simulations. Red lines indicate values from real cell tracks. Data points in left and right panels resutl from the same simulations.. . . 28
List of Figures vi
2.5 Effect of model parameters of preferred direction model on average speed and persistence time of simulated cells. In the preferred direction model the λDir mainly determines average speed of simulated cells and Pprst
mainly influences persistence times. Data presented are results from first and second grid search in single cell simulations. . . 29
2.6 Effect of model parameters of act direction model on average speed and persistence time of simulated cells. In the act model λDir and Aact both
influence average speed and persistence times of simulated cells. Data presented are results from first and second grid search in single cell sim-ulations. . . 29
2.7 Filling empty environment with non-moving cells. Green line represents theoretical limit of maximum total volume : 643 = 262144, and maximum number of seeded cells : 262144/150 = 1745.63 . . . 30
2.8 Effect of density on cell migration behaviour. A density increases in simu-lations of preferred direction model, cells move with higher average speed, higher persistence and in a coherent manner. Simulated cells in the act model move with lower average cells speed in denser environments, but also move more persistently and more coherent. . . 31
2.9 Results of parameter search of preferred direction model for full environ-ment including FRC network. Red line indicates goal values as observed in real cell track data. Average speed is accurately replicated by the pre-ferred direction model, but persistence times are consistently lower then real cell track data. . . 32
2.10 Fitting act model parameters to real data in fully populated environment including FRC structure. Persistence times in the act model are lower then in real cell track data, however, slightly higher then in the preferred direction model. Red line indicates goal values as observed in real cell track data. . . 32
2.11 Effect of graph theoretical topology of FRC structure on individual mi-gration behaviour of cells. . . 33
2.12 Turning angles of cells in preferred direction and act model in simulations with different FRC structures.. . . 34
2.13 Global order measure for collective movement in simulations with different FRC structures.. . . 35
2.14 Emergence of collective migration over time in simulations with different FRC structures. When no FRC structure is present all cells form one coherent stream. FRC structures with lattice-like properties allow for more coherent movement to form then FRC structures of more random topology.. . . 35
2.15 Relation between ω small-world network of FRC structures and measures for T cell migration behaviour in preferred direction model. . . 35
2.16 Relation between ω small-world network of FRC structures and measures for T cell migration behaviour in act direction model. . . 36
2.18 Comparing speed and coherent movement of cells in large gaps between simulated network to cells near centers of longest FRC fibers. Cells within a distance of 5 mum of gap centers and cells within a distance of 5.74 µm were sampled. Figure A shows the average number of cells within these distances. Figure B shows that cells move faster near in gaps then along fibers. However cells along fibers move more coherently in one direction, as represented by the increase order value in C. . . 37
2.19 Cell movement near centers of short fibers was observed to be more co-herent then cell movement along long fibers. Cells moving in larger gaps also moved slightly more coherent then cells in smaller gaps, however dif-ferences are very small when compared with coherent movement along fibers. . . 38
Abbreviations
LN Lymph Node
DC Dendritic Cell
LEC Lymphatic Endothelial Cell
BEC Blood Endothelial Cell
FRC Fibroblastic Reticual Cell
CPM Cellular Potts Model
ER-graph Erdos Renyi graph BA-graph Barabasi Albert graph WS-graph Watts Strogatz graph
GM-graph random GeoMetric graph
APC Antigen Presenting Cell
HEV High Endothelial Venule
ECM Extra Cellular Matrix
MCS Monte Carlo Step
0.1
The adaptive immune response and lymph nodes
The immune system is an immensely complex system that is of vital importance for human survival. The global pandemic of the SARS-CoV-2 virus is a strong reminder of the importance of an effective adaptive immune response. Lymph nodes (LNs) provide the environment in which many important processes take place that are vital for the adaptive immune response [2] [3].
T and B cells are key players in the adaptive immune response of the human body. B cells indirectly target foreign or malicious cells by secreting antibodies, proteins that bind specifically to their target cells, labeling them for attack by other cells in the immune system. Specific receptors on B cells bind to antigen, pieces of foreign organisms such as viral proteins, that cause the B cell to proliferate, differentiate into plasma cells and secrete antibodies. This differentiation of B cells into antibody secreting plasma cells often happens in germinal centers, specific regions of the LN that are mainly populated by B cells.
B cells also recruit helper T cells (CD4+) that play an import role in signaling via the secretion of cytokines. Some functions of this signaling include the stimulation of cytotoxic T cells and the differentiation of B cells into antibody secreting plasma cells. Another way by which the adaptive immune response can be initiated is via dendritic cells (DCs). They bind antigen at inflamed sites in the body and travel through the lymphatic systems until they encounter a helper T cell with cognate receptors to the specific antigen they are presenting. The proportion of T cells that have a receptor that can bind a specific antigen represented by a DC is about one in 105− 106 [4]. Hence it is very important that DCs encounter many T cells in a short period of time.
These encounters between DCs and T cells often happen in LNs, as a large amount of the T cell population resides there. Lymphatic fluid, containing DCs, B cells, T cells and other small blood molecules, enters the LN through afferent lymph vessels and flows into the subcapsular sinuses and to the LN cortex [5].The LN cortex is mainly populated by B cells, residing there until they get activated and form germinal centers. T cells mainly reside in the LN para-cortex, often called T cell zone, where they migrate in search of specific antigen presenting DCs [6].
Another type of cell that has more recently been found to be important for an adequate immune response is the stromal cell type. Different regions of the LN are populated by different subsets of stromal cell types. Stromal cells provide rigidity to the LN due to the collagen present in these cells. Stromal cells construct the sinuses, through lymphatic endothelial cells (LECs) and the vascular network, through blood endothelial cells (BECs). Stromal cells are involved in germinal center formation and B cell follicle
2
homeostasis in B-cell areas in the LN [7]. In the T cell zone the most abundant stromal cell type is the fibroblastic reticular cell (FRC), which constructs a network that guides T cell migration [8–10]. More specifically, T cells have been observed moving along fibers of the FRC network. Therefor the structure of the network is thought to provide a sort off ”road map” for T cell migration. Both T cells and DCs seem to stick to the FRC network, so FRCs are thought to mediate interaction between DCs and T cells [8]. Additional to directly influencing T-cell migration, FRCs also form a conduit system that transports small molecules such as antigen and signaling molecules such as chemokines and cytokines. Chemotactic gradients have been shown to influence T cell migration [11]. FRC cells have are known to express different types of ligands and receptors involved in immune cell functioning and homeostasis [12].
Although the above discussion shows the great advances being made in molecular im-munology, some questions remain hard to answer in vivo or vitro experiments. For exam-ple, cell staining used in experiments only stains a fraction of lymphocytes, which makes quantifying collective migration behaviour difficult. A more obvious disadvantage of wet lab experiments is the challenge of manipulating the structure of the FRC network. Novkovic and colleagues show that only after ablation 70% of the FRC, lymphocyte motility is altered [1]. However effectively changing the underlying graph theoretical topology of the FRC network has not been done yet. Also, because of the methodologi-cal challenges of recording T cell migration in the dense environment of the LN, recorded cell tracks are often relatively short, and the space sampled is limited.
Computational modelling of T cell migration can help in gaining more insight in the complex dynamics of T cell migration in the LN. In a computational model of T cell migration the virtual LN environment is easily manipulated and cell tracks of all simu-lated cells can be recorded and analysed. The Cellular Potts Model (CPM) has proven useful in simulating migration behaviour of biological cells [13,14]. The CPM allows for a wide range of movement patterns, combined with realistic morphology of simulated cells [15]. Cell shape of moving T cells is very flexible, constantly changing between very elongated and flat shapes, to more spherical [16]. This flexibility allows T cells to move through dens environments such as the LN or skin epithelia.
The migration behaviour of naive T-cells in the LN resembles brownian motion, with some variation in both speed and persistence of direction of movement. T cells have been shown to alternate periods of slow random motion with periods of faster more persistent motion [17–23]. This coupling between speed and persistence can accurately be modeled using CPM simulations [24]. T cells also show signs of collective migration by the formation of small streams of cells that migrate in same direction [25]. This stream formation behaviour has also been shown in simulations in the CPM [2]. It is
hypothesized that the formation of streams emerges because of the dense environment of the LN. T cells constantly ’push’ other T cells in random directions, until by change some groups of T cells align their movement and push neighbouring cells in the same direction, causing small streams of aligned migration to form.
The influence of the FRC network on this stream forming behaviour remains unknown, but as the FRC guides T cell behaviour it seems a good candidate in influencing stream forming behaviour in the LN. We suspect two possible ways in which the FRC network structure can influence T cell movement. As said earlier, T cell stick onto FRC fibers and migrate along them. This may cause streams of persistently moving T cells to form along side the fibers of the FRC network. A second possibility is that in the dense environment of the LN, FRC fibers obstruct the persistent movement of T cells, whereas in the gaps between fibers, T cells can move freely and might form coherent streams. The aim of current project is to inquire about the influence of the structure of the FRC network in CPM simulations on collective and individual migratory behaviour of T cells. To realistically model the structure of the FRC network in the CPM, some means of quantifying its structure is needed. Different measures from graph theory alongside some spatial measures will be experimentally applied for this purpose. This forms the subject of the first chapter in this report. Once an adequate model of the FRC network has been found, it will be used in CPM simulations to test the influence on T cell migration. One prerequisite for this is finding parameters in the CPM that allow for realistic simulation of T cell migration. To do so, cell track data from a intravital 2-photon experiment, profided by the lab led by Judith Mandl at McGill University in Canada, will be used to validate CPM parameters. This forms the subject of chapter two.
Chapter 1
Quantifying and 3D modeling of
the fibroblastic reticular cell
network structure
1.1
Introduction
The aim of this chapter is to find a robust way of generating a realistic 3D model of the FRC network. To do so, measures are needed that can be used to validate the generated models. Previously, graph theory has been used to explain the resilience to damage of the FRC network [1]. Moreover, the finding that lymphocytes move along the FRC network make the topology of the network a possible predictor for lymphocyte mobility. Graph theory has proven an useful means in many different domains to reveal intrinsic properties in many different application domains ranging from social networks [26], power grids [27] and biological networks [28,29].
Here we aim to use graph theory for two purposes. First we use different graph theoretical models as a blueprint for structurally different 3D models of the FRC network. Novkovic and colleagues use the notion of small-worldness to explain how the functionality of the FRC network is robust to damage [1]. A small-world network is characterized by short path lengths between any pair of nodes, yet also has relatively high clustering coefficients, and was introduced in the well known paper of Watts & Strogatz [30]. For T cell motility this would imply that T cells can migrate easily between any two nodes. However this hypothesis shows one shortcoming of using graph theory in this context; graph theory generally does not take into account positioning of nodes in space. Hence, graph theoretical and real life path lengths can greatly differ. To test weather small-world properties of the FRC network we aim to create FRC models varying in small-wolrd
properties. If the topology of the networks really is important for lymphocyte migration behaviour, then small-world measures should be able to predict certain aspects of T cell migration such as persistence of movement. Secondly, following the line of work of Novkovic and colleagues we aim to verify how well the generated 3D models resemble the real FRC network.
Although this notion of small-worldness is often mentioned qualitatively, it it can be quantified by comparing a given graph to an equivalent random graph. The small world measure σ is defined as follows :
σ = C/Crand
L/Lrand (1.1)
Where C, Crand, L and Lrand denote the average clustering coefficient and path length of the graph and of a random graph of with the same density and same number of nodes. Hence when a graph has a higher degree of clustering then a random graph and shorter path lengths, its σ value increases. σ Values higher then one indicate that a network is a small-world network [1].
One downside of this measure is that networks with very small clustering coefficient can easily be characterized as small-world, as long as they also have small path lengths. To counter this phenomenon another measure exists that compares a graph at hand not only to a random graph, but also to a regular lattice [31]. This measure ω is defined as:
ω = Lrand
L −
C
Clatt (1.2)
Networks of ω values near zero are considered small-world networks. Positive values of ω indicate more random properties in a network, and negative values of ω indicate a networks is more lattice-like [31,32]. The group of Novkovic reports a sigma of 6.7 and an omega of -0.27 for the graph based on the FRC network. Another graph theoretical measure that will be used to validate 3D models of the FRC network is the degree distribution of the networks.
The small-world property of the FRC network can likely be explained by geographical preferential attachment of fibroblastic cells during formation of the LN [33]. In their model Soekarjo and colleagues randomly seed simulated cells in a volume and allow for the formation of connections depending on a monotonically decreasing probability function. They show that this type of model greatly resembles the FRC network. This finding suggests another type of graph that can be used to generate realistic 3D mod-els of the FRC; a 3D random geometric graph. In a 3D random geometric graph, N nodes are randomly placed in a unit cube, with connections between nodes depending on a radius R. Instead of a monotonically decreasing probability function as used by
6
Soekarjo and colleagues,in a random geometric graph connections between nodes are formed depending on a Heaviside step function.
The exact mechanisms of how the FRC network is constructed during growth is still unknown. However Katakai and colleagues do propose a possible mechanism of FRC network formation [9]. They propose that the network formation is initiated by accu-mulation of lymphocytes and FRC progenitors, mediated by secretion of cytokines from lymphocytes. Subsequently, interactions between lymphocytes and FRC progenitors cause the FRC progenitor cells to generate fibers and form a connected network. The generation of 3D random geometric graphs resembles this proposed of FRC formation, and would therefor be a biologically plausible model.
The adaptive nature of FRCs cause the morphology of FRCs to vary depending on environmental factors. For example, during infection, the FRC network stretches out, giving rise to a dramatic expansion of the LN. This adaptive nature of FRCs also give rise to some variation between reports regarding the morphology of the cells. The morphology of a cell can be described through the morphology index which compares the surface area and perimeter of a given cell to circle as follows :
Imorph=
perimeter2
4π · area (1.3)
Hence, a perfectly circular cell would have a morphology index of 1. FRCs have been reported to have an average morphology index of around 8 [34], due to the long pro-trusions that are formed between cells in the network. Novkovic and colleagues use the notion of sphericity to describe the cell shape and report a sphericity of approximately .5 for FRCs [1]. Reports of cell area of FRCs range from 1 × 103 square microns to 2.7 × 103 square microns [1,34,35]. The average number of cell protrusions per FRC is reported to be between 4 and 5. The length of the protrusions of FRC range between 20 and microns +- 75 microns [1,34,35]. This variety in the shape of FRC cells contributes to the difficulty of generating a realistic model of the FRC network.
Above we mentioned the hypothesis that topology of the network determines T cell motility, based on in vivo experiments that show T cell moving along fibers of the FRC network. However, what these experiments do not show is that in fact the LN environment is completely filled with lymphocytes. The idea of singular T cells following a ’road map’ might depict the wrong picture. More accurate might be that a ’crowd’ of T cells, being pushed through the LN, in which the role of the FRC network would be more to mix the incoming T cells, and promote random search pattern. Following this analogy of crowd dynamics it is also likely that negative spaces, the gaps between fibers, determine T cell motility. To quantify the negative spaces in generated models
a gap analysis will be performed. We can compare the results of the 3D models to real data provided by the lab led by Judith Mandle at McGill University.
8
1.2
Experiments
1.2.1 Small-world properties in different graph models
Multiple graph models were used to generate graphs for the following purposes. First, we aim to find graph models that consequently generate graphs with small-world properties similar to that observed in the FRC network. Also we want to investigate how small-world properties vary in graph models that we expect not to be small-small-world. We used a wide variety of graphs including : Erdos-Renyi (ER) random graphs, Barabasi-Albert (BA) graphs, Watts-Strogatz (WS) graphs, random geometric graph (2D and 3D) and a variation on the Barabasi-Albert graph that allows for more tunable clustering (BA-cl) [36]. Sizes and densities are varied across simulated networks to see if found small-world values hold over larger and more dense networks. Graphs were generated of sizes 100,176,500 and 1000 nodes. Graph densities were varied by manipulating the expected average degree of the graphs. Expected average degrees (E[hki]) of simulated networks were : 16, 10.66, 8, 6.4, 5.2. Generating graphs with specific average degree requires finding some correct parameters for the specific graph type. For ER graphs with N nodes the p parameter resulting in average degree E[hki] is given by : 2d(N-1). For BA graphs the expected average degree is directly set by parameter M for hki = 2M . In WS graphs the expected average degree is 2 × avd. To find the correct K parameter for the BA-cl graphs the following equation needs to be solved : K2+ KE/d + N d = 0. Finally calculating the expected average degree for random geometric graphs is quite complicated analytically, and therefore a polynomial fit was used to find a relation between number of nodes and the number of edges.
For every parameter combination 10 graphs were generated. Computing ω and σ mea-sures requires generating random graphs and lattices with the same size and density. To minimize computing time a database was created with 100 random graphs and lattices for every graph size and density.
1.2.2 From graphs to 3D models
One downside of using graph theory to quantify and model the FRC networks is that most graph models do not contain any spatial information or measures; the positions of nodes are often disregarded. In order to use the generated graphs for modeling the FRC network, they need to be placed in 3-dimensional space, with the exception of random geometric graphs. This was done using a force directed graph drawing algorithm developed by Fruchterman and Reingold [37], which aims to place nodes in a space in such a way that edge lengths are evenly distributed. In this algorithm nodes are modeled
as particles, with repellent forces between two nodes, but attractive forces only between nodes connected by an edge [37]. Forces are calculated every time step using following equations : F a(d) = d2/k, F r(d) = −k2/d. k Here is the optimal distance difined as : k = C ∗ (area/N )1/3. The maximum displacement per time step is regulated by the
temperature, which is lowered every time step. Initial positions of nodes were random coordinates.
The Cellular Potts model, that is used for simulations in chapter 2, is a on-lattice model. Therefor, in order to use the generated spatial graphs, the positions of the nodes need to be discretized and edges need to be drawn in 3D space. We aimed for simulating a FRC network in a cube of 256 cubic microns, where approximately 17% of the voxels are occupied by FRCs, following the work of Beltmann and colleagues [4]. In all models one voxel corresponds to a cube of 1µm The thickness of connections and the number of nodes were varied in order to achieve this 17% occupancy. The random topology of some of the network types (ER,BA and BA-cl) cause the graph drawing algorithm to draw all nodes in a sphere. Subsequently placing such spherical FRC structure in a cube results in empty spaces in the corners of the cube, and highly cluttered center, which is would be a unrealistic model for the FRC network. Hence these graphs were drawn in a much larger space , 10243 microns, and a cube of 256 microns was cut out. Again the number of nodes of the graph was tuned to achieve 17% occupancy.
1.2.3 Gap Analysis
As an additional spatial quantification of the modelled FRC structures a gap analysis was performed. Random slices of 256x256x4 of the 3D models were taken and compressed into a image, similar to the procedure described by Acton and colleagues [38]. This picture was enlarged to a dimension of 1080x1080 to allow for a wider range of radii to fill the gaps. A euclidean distance transform was used to determine the largest gap between simulated FRCs, and the space was iteratively filled by circles with diminishing radii. For every radius size, the frequency was recorded and used to compute the proportion of space filled by a circles with a specific radius. In this way a relation between radius and the proportion of total surface area can be compared between models. Using data from an analysis performed by Shabaz Sultan in collaboration with Connie Shen from McGill University, the gap analysis of our FRC models was compared with gap analysis of pictures of FRC structures in mouse.
10
1.3
Results
1.3.1 Small world properties of different graph models
The aim of the first experiment was to find graph topologies with similar small-world properties as the FRC network, and then to assess whether small-world measures were consistent when network size and densities are altered. Across all types of networks created, random geometric and Watts-Strogatz graphs most frequently have similar σ and ω values to the FRC network as reported by Novkovic and colleagues. Figure 1.1
shows sigma and omega values of all generated graphs, where individual dots represent different combinations of size and average degree. For the Watts-Strogatz and BA-cl graphs additional tunable parameters cause a wider range of sigma and omega measures.
Figure 1.1: Small-world measures for all types of simulated networks. Observed small-world properties of FRC as observed by Novkovic and colleagues were σ : 6.7 and ω : −0.27 .Points represent average σ and ω values for generated graphs of certain size and density. In the WS and BA-cl models one extra parameter can be varied, which causes the higher amount of data points in these graphs. Size of the points show summed variance within graphs of specific density and size, ranging between .0054 and
1.345.
Graphs that most closely resembled the σ and ω value of the FRC were : random geometric graphs (N = 176, hki = 8, r = 0.27, σ = 7.76, ω = −0.26), Watts-Strogatz
graphs (N = 176, hki = 8, P = 0.077, σ = 6.6, ω = −0.24) and BA-cl graphs (N = 100, hki = 6.4, σ = 6.6, ω = −0.31)
However we must note that the method used to calculate σ and ω values was not completely accurate, see table1.1. This could either be due to stochastic nature of some of the graph models, or due to mistakes in calculating the expected average degree. As discrepancies due to stochasticity should cancel out over multiple graphs, we conclude some mistakes must have been made in calculating parameters of different graph models. As a result the average degree of generated graph models differs from the random graphs and lattices that they are compared to in order to calculate σ and ω. The exact effects there of however have not been tested.
Table 1.1: Percentage of absolute error of average degree per graph type. As graphs are compared with random graphs of same size and density to compute ω and σ, error
in estimated average degree causes error in these measures.
ER BA BA-cl WS GM-2D GM-3D
0.36% 11.59% 11.92% 9.33% 10.82% 32.76%
(a) σ – network size (b) ω – network size
(c) σ – network density (d) ω – network density
Figure 1.2: Effect of networks size and density on σ and ω small-world measures. Sizes of the dots indicate standard deviation of σ and ω measures. σ Small-world measure is clearly influenced by both network size and density. Whereas ω measure captures differences between networks more consistently across different network sizes
and densities.
Figure 1.3 shows the influence of network size and average degree on both small-world measures. In general the omega value seems less influenced by network size and density,
12
and there is a more consistent separation between network types in omega values. For most graph types sigma values increase as network size grows, and decrease as the average degree of the networks increases.
The degree distribution of the FRC network has been reported to be normally distributed with an average of 7.78 [33]. Figure 1.3shows the degree distributions for graphs with an expected average degree of 8. Figure 3 shows that this distribution is most closely matched by a Erdos-Renyi network. However the degree distribution of 3D random geometric graph also seems a close fit, except for the positively skewed shape, see figure
1.6.
(a) Erdos-Renyi (b) Barabasi-Albert
(c) BA-cl
(d) Wats-Strogatz (e) 3D Geometric
Figure 1.3: Degree distributions for generated graphs with similar size and density as FRC network reported in [1] : N = 176 , hki = 8. Specific parameters of graph models : PER = 0.0457, MBA = MBA−cl = 4,PBA−cl = 0.666, KW S = 8,PW S = 0.0278,
rGM −3D= 0.1234.
Edge lengths of the FRC network have been reported to be follow a positively skewed distribution around an average length of +- 30 microns [1,34,35]. Figure 1.4shows the
edge length distribution of the simulated networks. Erdos-Renyi, Barabasi-Albert, Wats-Strogatz and BA-cl networks show a high amount of long range connections compared to the FRC network. The edge length distribution of 3D geometric model is negatively skewed in contrast with the FRC network. Overall the edge distribution of the WS networks seems closest to the FRC as it is positively skewed with a mean length of 42.1 microns.
Figure 1.4: Distribution of edge lengths of 3D models based on different graph types.
Visually however the GM networks seems most resemble the structure of the FRC net-work best, which is nicely illustrated in figure1.5. On the left a cross section of the 3D random geometric model, in the middle an actual cross section of the FRC network. In
1.5b the yellow and light green lines represent fibers of the simulated FRC structure. The remaining space is filled with circles, which is an example of gap analysis used in attempt to quantify the visual resemblance.
(a) Real rat LN (b) Simulated FRC
Figure 1.5: Cross-section of stromal structure in mouse LN and simulated stromal cell structure, filled with circles after gap analysis. The yellow structure in B represents
the simulated stromal structure.
However results of the gap analysis do not reflect the visual resemblance between sim-ulated structure and real networks. In figure 1.6 shows the proportion of surface area filled by circles of a certain radius. In LNs of real mouse, larger gaps are far more present
14
then in any of the graph based models. In mouse LNs gaps are present that can be filled with circles with radii up to 56 µm, whereas the largest radius observed in any of the graph based models was 17.25 µm in the GM model. Figure 1.6 also shows that the spacing between FRC fibers can vary with age of mouses. To find out which model best fitted the mouse data, the mouse data was averaged over young and old mouse. Then the error was computed of interpolated curves of the graph models. Interpolation was done by using B-spline representation of the raw curve over the same interval, with a degree of 3 and smoothing condition : P
Y(y − g(x))2 <= s,where g(x) is the smoothed interpolation, Y the set of data points and s was set to 100. Error of real data points could not be computed as data points did not align. The random geometric model had lowest percentage absolute error, see table 1.2
Figure 1.6: Proportion of surface area filled with circles of varying radii in real mouse FRC networks and graph based 3D models. Gray lines represent variety of gaps present in real mouse LNs. Left pannel shows raw data, right pannel smoothed data using B-spline representation. Error between smoothed curves of grahp models and mouse data
was lowest for the random geometric model (GM-3D).
Table 1.2: Percentage of absolute error between interpolated curves from models and real mouse data.
ER BA BA-cl WS GM-3D
1.3.2 Conclusion
The small-world properties of the FRC network seem most consistently replicated by random geometric graphs. However as densities of graphs were not correctly estimated, especially for the random geometric graph, this result is not reliable. Out of both small-world measures computed, ω seems to capture structural differences between graphs best, independent of graph size and density. If small-worldness of the FRC network does influence the migration behaviour of T cells, we expect the ω value of a network to be able to predict some properties of T cell movement in simulations.
Although the degree distribution of the random geometric graph resembles the real FRC reasonably well, long range connections that are present in real FRC networks are not modeled accurately by the random geometric graph. A more accurate model would include more long range connections similar to the model described by Soekarjo and colleagues [33]. However based on the results of the gap analysis and the notion that it is the most biologically plausible model, the random geometric graph will be used as model for the FRC network in the next chapter.
Chapter 2
Simulating migration behaviour
of T-cells in lymph node
2.1
Introduction
In this chapter we aim to realistically model T cell migration behaviour in a densely populated LN, and assess the role of stromal cell structure on collective and individual migration behaviour.
The main objective of T cell migration within the LN is searching for antigen pre-senting cells (APCs). The most effective search strategy for finding APCs would be to move directly towards them. Antigen presenting cells enter at specific sites (HEVs and sinuses), hence if all lymphocytes would move directly towards them general lymph flow would be greatly impaired. Also, only one T cell in 105 - 106 is able to bind to a specific antigen [39–41], and it is therefore of great importance that many T cells make contact with an APC. Observations suggest that T cells accomplish these goals by switching between short periods of high speed movement into a consistent direction and slow movement with sudden random turns. On longer time scales T cell movement is seemingly random and this movement pattern has been described as ”stop and go” movement [17–20, 22]. Periods of persistent movement allow cells to cover distance on longer time scales, enabling them to circulate through the body. The periods of more slow and random movement allow T cells to scan areas in search of APCs.
The FRC network is thought to facilitate the ”stop and go” behaviour as migrating T cells are seen moving along its fibers. This suggests the that the FRC networks acts as a ”road-map” for T cell migration, making the topology of the FRC network a possible factor in determining T cell motility behaviour. However, simulations suggest
that the density of the LN environment influences motility of simulated T cells [4]. We hypothesize that in this dense environment it is actually the gaps between FRC fibers that allow for streams of fast and persistently moving T cells. Beltman and colleagues use randomly placed rods as FRC fibers in their experiment, filling 17% of simulated space. As shown in the gap analysis of previous chapter, gaps between fibers of more random models are relatively small (see figure 2.1a as example of more random FRC structure). We suggest that structures with larger gaps, more lattice like structures, allow for more persistent movement of T cells, and more coherence between moving T cells in these gaps (see figure2.1b as examples of lattice like structures. It is important to note that lattice like networks have ω values close to minus one. Gaps in lattice like structures will also be more uniformly distributed. In this sense the gap analysis and small-world measure ω can capture the same property of a network or structure.
(a) Erdos-Renyi based FRC structure (b) FRC structre based on a random geometric graph.
Figure 2.1: Example of more random like FRC structure (A) and more lattice like FRC structure (B).
Due to the random properties of T cell movement, it is often quantified by comparing it to Brownian motion. One method often used is plotting the mean square displacement of a cell track as a function of time. For Brownian motion this relation is linear. From this relation a diffusion coefficient and persistence time can be estimated by fitting the curve to F¨urth’s formula [42]:
h ~d(t)2i = 2nD(t − P (1 − e− t/P)). (2.1)
Here t is time, n the dimension of the space in which a cell moves, D the diffusion coefficient, often called motility coefficient for cells, and P the persistence time. This estimation of diffusion coefficient and persistence time through this formula however introduces extra variation caused by the fitting procedure.
18
Another way of quantifying how persistent a cell moves is by computing autocorrelation between instantaneous movement vectors at different time intervals within the same cell track [24]. This offers a more intuitive way of quantifying the persistence of a cell track, as the correlation between two vectors at different time points corresponds to the cosine of the angle between the vectors. When vectors are pointing in the exact same direction the cosine between them is 1, and vectors opposing each other have cosine of -1. A completely random cell track will therefor have a average autocorrelation of 0 for any time interval and a cell moving straight in one direction will have a autocorrelation of 1 for any time interval.
In the dense environment of the LN, T cells are packed together in a crowd-like manner, yet remain highly motile. The FRC network has been shown to influence movement of T cells in two photon microscopy experiments in vivo [8]. However these experiments can only record a small proportion of the population of T cells present in the LN, and therefore do not show the dynamics of T cells moving as crowd. Figure 2.2 illustrates this difference in a simulated environment. In this chapter we aim to quantify crowd dynamics and the influence of the FRC network in a simulated environment.
(a) Subset of simulated cells in CPM sim-ulation.
(b) Example of CPM fully populated by simulated T cells.
Figure 2.2: Difference between only observing a subset of cells, versus the fully popu-lated environment of the LN. Thin lines in the left panel show represent the track of the center of mass of the simulated cells. Gray structure is the a model of FRC structure based on random geometric graph. Different colors of cells have no significant meaning.
Dimension of simulated environment is 32 cubic microns.
The most straightforward way of quantifying whether cells move coherently in a similar direction is summing the normalized instantaneous velocity vectors at a specific time of all simulated cells. When movement is random, the sum of these vectors would be zero, and the sum grows as more cells move in similar direction at a specific point in time. When all cells move in a similar direction this sum should approach one. Turning angles can also be useful to illustrate variation in cell movement, showing how much turbulent
the crowd of simulated cells moves. Large average turning angles in a simulation indicate that cells move in a more Brownian manner, and variance in turning angles can show how big differences in persistence are between cells in a simulation.
To model T cell migration, two variants on the Cellular Potts Model (CPM) will be used. The CPM is a on-lattice model in which a cell is modelled as a group of particles with the same ”spin” (the exact model definitions are described in the following section). Movement of the cells is modelled using either the preferred direction model [4] or the act model [15]. In the preferred direction model, movement is forced upon the cell, whereas in the act model it emerges through a memory system inspired by actin remodelling in real cells. Niculescu and colleagues identify two modes of cell movement in the act model. For relatively low act values simulated cells resemble movement patterns of amoeboid cells, where wave fronts of high act values subsequently emerge at different sides of the cell, causing the cell to stop and move in random directions. For higher act values the cells form a consistent wave front at one side of the cell, causing the cell to flatten and move consistently along its longest axis. They categorize this second type of movement as keratocyte like. In the preferred direction model, the persistence of cell movement can be directly set through model parameters whereas in the act model it results from memory of lattice sites of cells. The cell movement in the act model is therefore a lot harder to predict, and realistically modeling cell motility in a densely populated 3D environment has not yet been done. In the preferred direction model it is more likely to find a good fit to real data, whereas the underlying mechanism that allows the cell to move is more realistic in the act model. Preferably we would use the act model, but as it is not certain yet that parameters can be found that result in realistic cell motility in such a dense environment, we will use both models.
The main reason the CPM is used in this project is because the morphological flexibility that can be modelled in cells using the CPM. Especially for T cell the flexibility of their shape plays an important roll in motility, as it enables them to crawl through small gaps [16]. Other models used for random motion of particles do not include this flexibility of shape.
To validate realistic movement of simulated T cells, model parameters will be fit to real cell track data from the the group of Judith Mandl and colleagues [43], made publicly available trough Motility Lab [44]. This data contains cell tracks of both CD4+ and CD8+ T cells traversing through mouse LNs. Before we can start simulations with dif-ferent stromal structures, we first need to establish that our model can replicate realistic T cell movement. Hence, as a first step, single cell simulations in an empty environment will be used to ensure both models can simulate realistic cell behaviour. Once good parameters are found, simulations with increasing densely populated environments are
20
performed to find how the density of simulations influences the motility behaviour of simulated cells. Finally we will search for realistic cell behaviour in a fully populated environment including stromal cell structure. Then by manipulating the structure of stromal cell network, the influence of the network on migration behaviour of T cells is researched.
2.2
Methods
2.2.1 Cellular Potts Model
The cellular Potts model (CPM) is a on-lattice model in which cells are comprised of several grid cites of same cell type (τ ) and identification number (σ). The simulation advances by means of random copy attempts of a certain cell to neighboring grid cells. The success of the copying attempt is determined by the effect on the global energy of the system; which is described by the Hamiltonian:
H =X u,v Jτ (σu),τ (σv)(1 − δσu,σv) + X σ λArea(aσ− Aσ)2+X σ λP erimeter(pσ− Pσ)2. (2.2)
The hamiltonian consists of three terms. The first part describes the interaction energies between neighboring grid cells (u and v) of different cell IDs, using the Kronecker delta function (δσu,σv) which is one when σu = σv and zero otherwise. The second term
makes sure that the total volume (aσ) of a cell stays close to the goal volume of that cell (Aσ). The last term has a similar function as it controls the possible cell shapes and ensures that the collection of grid cells compromising a single cell stays connected. When a copy attempt decreases the overall energy of the system it is always accepted, otherwise the chance of accepting the attempt is equal to the Boltzmann probability e(dH/T ), with ∆H the difference in energy and T the temperature of the system. A Monte Carlo step is defined as the number of copy attempts equal to the number of lattice sites. Different cell types are created by defining setting different parameters for different types. Energies between cell types determine the how much cells are attracted or repulsed to one another, and differences in cell area and perimeter parameters allows for heterogeneity in morphology between cell types. In simulations we use a minimum of two cell types, the background or extra cellular matrix (ECM) and T cells. To make the T cells motile two different models are used; the preferred direction model and the act model. All models were implemented using C and Python programming languages.
2.2.2 Preferred direction in CPM
In the preferred direction model every cell has a individual preferred direction vector ( ~P ) and a persist (Pprst) parameter that is set for all cells of its type. The persist parameter determines the amount of directional diffusion to ~P every MCS following equation:
~
22
Where ~Pt is the current preferred direction vector of a cell, ~Pt−1 is the normalized direction vector of a cell 15 MCS ago, and ~Vt is the normalized displacement vector between the center of mass of a cell at this MCS and the previous MCS. When calculating the energy difference (∆H) of a copy attempt of a certain cell, a term is added that gives copy attempts in the same direction as ~P a higher probability. This term is described by :
∆Hv = −λDir× ( ~P · ~XC→v). (2.4)
Here λDir ss the λ energy parameter set for that cell type and ~XC→v the normalized difference vector between lattice cite v and the current centre of mass location of the cell.
2.2.3 Act model in CPM
The act model for cell movement depends on a memory mechanism and random fluctu-ations in cell membrane present in CPM simulfluctu-ations. Movement in this model emerges when over a period of time more protrusions happen at a particular side of a cell, causing a bias for cell movement in that direction, resembling actin remodeling that drive cell movement in real cells. When a given cell of protrudes into ECM, that newly occupied grid cell is given a value; the maximum act value (Aact) which is set as a parameter. Every MCS this value diminishes by one. The collection of act values present in a cell biases protrusion close to sites with high act values. This is done by subtracting the following term from the ∆H calculation :
∆H = λDir
(GMact(u) − GMact(v)) (2.5)
GMact(u) is the geometric mean of act values for lattice site u, calculated through : GMact(u) = ( Y
y∈V (u)
Aact(y))1/|V (u)|
(2.6)
V (u) denotes the direct Moore neighborhood of grid cell u and Aact(y) the act value at lattice site y. In the act model cell movement occurs when a front of high act values form inside a cell, biasing copy attempts in that direction.
By varying parameters Aactand λDira wide range from persistent to random cell move-ment can be modelled [15], with realistic cell morphology, emergent from a local intra-cellular mechanism.
2.3
Experiments
A set of different experiments were conducted with multiple aims. First single cell simulations were performed to assess the relation between simulation parameters and cell motility behaviour and to confirm that both models can be used to model real-istic individual T cell motility behaviour. Then experiments were performed with an increasing number of simulated cells to see how the behaviour changes as the simulation environment becomes more densely populated. Subsequently we tried to find model pa-rameters resulting in realistic cell motility in a fully populated environment, including a FRC structure based on a random geometric graph. Then simulations with varying FRC structures were performed to see how the FRC structure influences the motility of simulated cells. In all experiments periodic boundaries where used.
2.3.1 Single cell simulations
To find model parameters resulting in realistic cell movement, simulations of a single cell surrounded by ECM were performed. Cells with a target area of 150 voxels were seeded in a cube of ECM of 64x64x64 voxels. The size of one voxel corresponded to 1 µm, resulting in cells with a diameter of approximately 3.3 µm. The λDir parameter of both movement models was varied from 50 to 20.000, with regular intervals of 100 between values 100 and 2000. For the act model,Aact values were varied between 10 and 250 in 25 regularly spaced steps, the Pprst parameter of the preferred direction model was varied from .1 to 1.0. Every parameter combination was repeated 10 times. For overview of other CPM parameters see table2.1. Most of these parameters were chosen based on previous work [4,24]
Results of this first grid search were compared with real data to find a range of param-eters that provided a starting point for the second more fine grained grid search. For the act model λDir parameter was varied from 100 to 2000 in 20 steps and the Aact pa-rameter was varied between 40 and 80 in 21 steps. The Pprst parameter of the preferred direction model for the second grid search ranged from .75 to .95 in 21 steps and the λDir energies ranged from 2000 to 4000 in 21 steps. After the second grid search the absolute error of every parameter set was evaluated to find the best fit for both models.
2.3.2 Simulations at different densities
To assess the effect of the density of LN environment on migration behaviour we sim-ulated cells in increasingly dense popsim-ulated environments. Both the act and preferred
24
direction model were used, with parameters that were found to best fit real cell be-haviour in experiment one. All other parameters were consistent with experiment one. One adaptation being that the grid size in current experiment was smaller than before: 32x32x32. Densities varied from .1 to 1.0 with an increment of .1, were a density of .1 consisted of 21 cells being seeded and a density of 1. consisted of 218 cells filling all available space in the simulation.
2.3.3 Simulating fully populated environment
To research the main question, how stromal structure influences T cell migration be-haviour, we first searched for realistic cell behaviour in a fully populated environment including stromal cell structure. The adhesion energy between simulated FRC fibers and T cells was set to -5, to mimic the behaviour of real T cells sticking to the FRC network. The stromal cell structure was modelled using a random geometric graph with 170 nodes, placed in a cube of 64x64x64 with connections between nodes up to 20 mi-crons apart, similar to procedures used in chapter one. The stromal structure populated 17.8 percent of the total space. The remainder of the space was populated by 1362 cells. To find good parameter values for λDir energies, Pprstand Aact parameters, a grid search was performed, similar to that of single cell simulations described above. Table
2.1shows the other CPM model parameters.
Table 2.1: CPM parameters in different simulations of parameter search for single cells simulations and fully populated environments, in the preferred direction model
and act model.
Parameter Single preferred dir. Single act Full preferred dir Full act
Temperature 200 200 200 200 Sampling every 50 MCS 40 MCS 100 MCS 40 MCS Simulation time in MCS 104 8 ∗ 103 2 ∗ 104 8 ∗ 103 Burn in 100 MCS 100 MCS 50 MCS 50 MCS Volume (voxels) 150 150 150 150 λV olume 25 25 25 25 Perimeter 1200 1200 1200 1200 λP erimeter .2 .2 .2 .2 Adhesion T-ECM 0 0 0 0 Adhesion T-T - - 10 10 Adhesion T-FRC - - -5 -5
2.3.4 Effect of different FRC structures
To investigate the effect of different stromal cell structures on individual and collective migration behaviour, migrating T cells were simulated in full environments including
stromal cell structures with different topologies. Both act and persistent direction models were used, with parameters found in above described grid search. For the act model these parameters where a lambda energy of 1200 and a max act value of 40. For the persistent direction a lambda energy of 600 and persist parameter of .9 were used. The stromal structures were created using the same method as in chapter one. The number of nodes of the graphs were chosen to ensure that all modeled structures occupied approximately 17% of the environment. See figure2.3for 3D representation of the used structures. For every FRC structure five simulations were run of
(a) Erdos-Renyi (b) Barabasi-Albert (c) BA-cl
(d) Wats-Strogatz (e) 3D Geometric
Figure 2.3: FRC models based on different graph types.
2.3.5 Measures to quantify migration behaviour
Individual cell migration behaviour was quantified using two measures: speed and per-sistence time. Speed is calculated by summing all norms of the instantaneous speed vectors of a specific cell track :
V = 1 N · T X n∈N X t∈T |~Vn,t| (2.7)
Where N is the number of cells in and T the total of sampling points in the cell track. The persistence of a cell track is calculated based on the autocorrelation of a cell track, using the same method described by Inge Wortel and colleagues [24]. For every time interval δt, the correlation between instantaneous speed vectors δt sampling points apart is calculated within every cell track, which corresponds to the cosine of the vector at t
26
and vector at δt, and is calculated as follows :
r = Vt~ · ~Vt+δt
||~Vt|| · ||~Vt+δt||. (2.8)
Then the time interval is estimated using a finite difference approximation at which the autocorrelation has decreased to a value of .5, i.e. persistence half time (Pht) :
Pht= δtr>.5+
.5 − rδtr>.5
(δtr<.5− δtr>.5) · (δtr<.5− δtr>.5). (2.9) Where δtr>.5 and δtr<.5 are the time intervals for which the correlation is closest to .5. If the autocorrelation for every is higher then 0.5 the persistence time is estimated to be the largest , which is set to 100. Persistence half times of 0.5 indicate that cells moving in a highly random manner, whereas persistence times of 100 indicate cells move consistently in one direction. Estimating persistence half time can be done individually per cell, or can be calculated by computing an autocorrelation curve for all cells in a simulation together. Cell tracks resulting from in vivo experiments of real lymph nodes often differ in length, and hence the latter approach is often used to get a more accurate estimate of the autocorrelation for every δt. Hence we will used this in analysis of both simulations and real data.
For collective cell behaviour we compute two measures to quantify the amount of coher-ence in the migration behaviour between cells. The first one is the global order (Oglob) measure, which is computed by summing the normalized instantaneous speed vectors every time step in the simulation, and taking the average of the norms of this summed vector at every time step. Note that if all cells move in the same direction, the global order is equal to the average speed.
Oglob = 1 N · T X t∈T ||X n∈N ~ Vn,t |~Vn,t| || (2.10)
The other measure we compute for coherent migration is local order, which is very similar to the global order only here the space is divided sub regions and only the vectors within certain region are summed.
In the last experiment, where we compare migration behaviour of simulated T cells in presence of different stromal structures, we aimed to quantify weather streams of cells form in big gaps between the network, or along the fibers. To do so we performed a gap analysis in 3D on the stromal structures, and computed the order of cells within a radius of 5 µm of the ten largest gaps. We also computed order of cells within a radius of 5.74 µm of the centers of the largest fibers. We used a larger radius near the fibers
to compensate for the space occupied by the fiber itself. Finally to asses whether the size of gaps or the length of fibers can predict how coherent cells move, we correlated the gap size and fibers length. For fibers length we now sampled randomly, to ensure enough variation in edge length.
28
2.4
Results
2.4.1 Single cell simulations
Average speed and persistence time calculated from in vivo two photon experiments performed by Judith Mandl and colleagues [43] were 5.10 µm/min and 6.15 minutes, respectively. For both movement models good fitting parameters were found in the second grid search. For every parameter combination the absolute error was calculated to find optimal parameters. For the preferred direction model best fitting parameters were λDir = 2000 and Pprst = .87, with an absolute error of = 2.50, resulting in an average speed of 4.81 µm/min and persistence time of 5.89 minutes. Best fitting act parameters were λAct= 1100 and Aact= 44 with an absolute error of = 2.14, resulting in average speed of 5.32 µm/min and persistence time of 6.01 minutes. Figure2.4shows the speed and persistence for every parameter combination of this grid search.
(a) Preferred Direction model
(b) Act model
Figure 2.4: Average speed and persistence times as resulted from second grid search for single cell simulations. Red lines indicate values from real cell tracks. Data points
As expected the influence of λDirand persist parameters in the preferred direction model show a direct influence on speed and persistence. The Pprstparameter has little effect on speed and λDir seems to only amplify the effect of Pprst on persistence time, see figure
2.5. The persistence time of simulated cells increases exponentially with Pprst. This is a direct result of the implementation of the preferred direction model, as movement is forced on the simulated cells, and is less realistic then the act model.
Figure 2.5: Effect of model parameters of preferred direction model on average speed and persistence time of simulated cells. In the preferred direction model the λDir
mainly determines average speed of simulated cells and Pprst mainly influences
persis-tence times. Data presented are results from first and second grid search in single cell simulations.
Figure 2.6: Effect of model parameters of act direction model on average speed and persistence time of simulated cells. In the act model λDir and Aact both influence
average speed and persistence times of simulated cells. Data presented are results from first and second grid search in single cell simulations.
As can been seen in figure 2.6, in the act model both the λDir and Aact influence the persistence time of the cell movement. The average speed of simulated cells in the act model is influenced most by Aact parameter instead of the λDir parameter. The Aact parameter also influences persistence times of simulated cells, as seen in the left panel
30
of figure 2.6. Previously the authors of [15] identified two types of movement possible in the act model, which could relate to the effect of Aact on persistence times as seen in
2.6on the right panel. Figure2.6shows that even though act model is biologically more relevant, model parameters influence cell motility in a non-trivial way. Presumably this makes finding parameters that result in realistic cell behaviour more difficult in simulations of more dense environments.
2.4.2 Increasing density of simulated environment
Before studying the effect of density on cell movement we tested whether non-moving cells would maintain size as the number of cells seeded in the environment increases. Cell sizes decrease slightly before the environment is theoretically full, see figure2.7.
Figure 2.7: Filling empty environment with non-moving cells. Green line represents theoretical limit of maximum total volume : 643 = 262144, and maximum number of
seeded cells : 262144/150 = 1745.63
When simulating moving cells in increasingly dense environments, the density does in-fluence cell motility. Surprisingly density has an opposite effect on speed in both models. As the simulation environment becomes more densely populated with cells, the average speed increases in the preferred direction model, whereas it decreases in the act model. In the act model cell movement is driven by an intrinsic memory mechanism based on fluctuations in cell shape. In dense environments these fluctuations might decrease and cells might obstruct each other in moving in a certain direction once a wave front of high act values has formed, possibly causing the decrease in average speed as seen in figure
2.8. The increase in average speed in the preferred is a quit surprising result, but might be attributed to the formation of streams as observed by Beltman and colleagues [4]. Possibly when cells move in streams the energy needed to move decreases, resulting in higher average speed.
Influence of density on persistence time is quite similar in both models, both show a large increase in persistence for densities above +- 0.3. At lower densities, interactions between simulated cells seem to not influence persistence times, whereas it does already
influence speed. Both models show a sharp, phase transition like increase in persistence times from goal value of around 6.15 to 100, which in this setup is the maximum value for persistence time. This phase transition can be seen as a transition from ’gas-like’ behaviour,wherein cells can collide but move freely most of the time, to ’liquid-like’ behaviour were cells continuously interact.
Figure 2.8: Effect of density on cell migration behaviour. A density increases in simulations of preferred direction model, cells move with higher average speed, higher persistence and in a coherent manner. Simulated cells in the act model move with lower average cells speed in denser environments, but also move more persistently and more
coherent.
Both global and local order measures increase with density in both models. Cells move slightly more coherent in the preferred direction model, again in line with previous work of Beltman [4]. The high amount of ordered movement present in these dense simulations is also likely to be caused by the relatively small simulation domain, and the periodic boundary conditions.
2.4.3 Simulating fully populated LN
Subsequently a parameter search was performed in a fully populated LN including FRC structure based on a random geometric graph. Finding model parameters that resulted in persistence times as observed in real cell track data proved difficult in both models. The majority of parameter combinations tested in the preferred direction model resulted in an average persistence time much lower then in real data. Both average speed and persistence times in simulations varied in a less predictable manner then in single cell simulations. Best fitting parameters in the preferred direction model were λdir = 600 and Pprst = .9 with an absolute error of = 9.92, resulting in a average speed of 4.75 µm/min (5.10m/min in emperical data) and persistence times of 0.90 minutes (6.15 minutes in emperical data).
Results of the parameter search in the act model were similar as most parameter combi-nations resulted in low persistence times. The relation between speed and both motility
32
Figure 2.9: Results of parameter search of preferred direction model for full environ-ment including FRC network. Red line indicates goal values as observed in real cell track data. Average speed is accurately replicated by the preferred direction model,
but persistence times are consistently lower then real cell track data.
parameters of the act model was similar as observed in single cell simulations. Param-eters best fitting to real cell track data were λact = 1200 and Aact = 40. with an error of = 9.26, resulting in an average speed of 4.48 µm/min and persistence time of 1.64 minutes.
Figure 2.10: Fitting act model parameters to real data in fully populated environment including FRC structure. Persistence times in the act model are lower then in real cell track data, however, slightly higher then in the preferred direction model. Red line
indicates goal values as observed in real cell track data.
2.4.4 Effect of FRC network topology
Even though persistence times in previous section were below that of experimental data, we were still interested in whether cell motility could be influenced by varying the FRC network. To do so, best fitting parameters from previous section were used in simulations
