Two mathematical approaches to a study of T cell motion and activation in the lymph node by Monica Delgado Carrillo B. Mathematics, Universidad de Guanajuato, 2009 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Mathematics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2012 © Monica Delgado Carrillo 2012 Abstract T cells are part of the immune system and as such play a very important role in keeping us healthy. One crucial step in the complex process which is the immune response to pathogens is T cell activation. The general goal of my thesis is to mathematically describe the migration patterns followed by T cells while waiting to be activated in the lymph node. Insight into these migration patterns could lead to better knowledge of the strategies T cells take to make activation such an efficient process. In order to fulfill my goal I have used two different approaches: one mainly computational and the other mainly theoretical. On the computational side, I analyzed three-dimensional microscopic movies of mice lymph nodes inside of which labelled T cells are moving. From the movies I extracted the trajectories of the cells. I studied movies from two experimental frameworks, exogenous and endogenous. On the former, more frequent type of experiment, T cells are labelled outside the mouse and then transferred in. The endogenous experiments, on the contrary, involve genetically modified mice whose T cells are born labelled. I concluded that there is a significant difference in labelled T cell motion between the two experimental frameworks. This suggests that previous results from exogenous experiments should be treated with caution due to possible errors introduced by the methods specific to that type of experiment. On the theoretical side I studied the time it takes for a model T cell to be activated under different scenarios regarding the characteristics of the lymph node as well as of the other cells in it. Since T cells become activated after establishing contact with a specific cell among many similar ones which also move within the lymph node, what I effectively computed was the mean first passage time for a model T cell to reach a defined target within the model lymph node. ii Preface My supervisor, Daniel Coombs has been involved in all areas of this thesis. Raibatak Das, Coombs’ postdoctoral researcher, was involved in the development of Chapter 2. They provided guidance and feedback, edited manuscripts, assisted with computational difficulties and gave advice in interpreting results. Chapter 2 is the result of work in collaboration with the laboratory of Dr. Pasquale Maffia from the Institute of Infection, Immunity and Inflammation at the University of Glasgow, in Scotland, UK. Dr. Maffia and his collaborators own the copy rights to the data. All experiments were performed by Ross McQueenie. McQueenie also helped in describing the experimental methods for Section 2.2.1. I did the analysis, aided by Das and Coombs. Some of the codes I used were taken from [13] either in full or in parts. I used MATLAB for all computational results from this chapter. The idea for the algorithm described in Section 2.4.3 was suggested to us by Khuloud Jaqaman. For Chapter 3 I received help and guidance from Professor Michael Ward. Figure 1.1 is here reproduced by kind permission of the Nature Publishing Group (License Agreement number 2970990947701). Figures 2.1 and 2.2 are here reproduced by kind permission of Michael W. Davidson of the National High Magnetic Field Laboratory, the author, under the condition that this thesis will only be available online in PDF format. All other figures were designed and created by me. iii Table of contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii List of tables Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Biological description of the problem . . . . . . . . . . . . . 1.2 Current approaches . . . . . . . . . . . . . . . . . . . . . . . 2 Analysis of imaging data from T cells . . 2.1 Two-photon microscopy . . . . . . . . 2.2 Experiments . . . . . . . . . . . . . . . 2.2.1 Experimental methods . . . . 2.3 Methods for the inference of cell paths 2.3.1 Single particle tracking . . . . 2.3.2 Steps in the tracking process . 2.3.3 Denoising with wavelets . . . 2.3.4 Track extraction . . . . . . . . 2.4 The detection step . . . . . . . . . . . 2.4.1 Spot detector . . . . . . . . . . 2.4.2 My detection technique . . . . 2.4.3 Using u-track for detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 2 4 4 6 7 7 7 9 10 10 14 14 15 17 iv Table of contents 2.5 2.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 22 25 31 34 First passage time of a T cell to an APC . . . . . . . . . 3.1 First passage time concepts . . . . . . . . . . . . . 3.2 Basic problem . . . . . . . . . . . . . . . . . . . . . 3.2.1 Neumann boundary . . . . . . . . . . . . . 3.2.2 Dirichlet boundary . . . . . . . . . . . . . . 3.2.3 Robin boundary . . . . . . . . . . . . . . . 3.3 Multiple targets with asymptotically null volume 3.3.1 Neumann boundary . . . . . . . . . . . . . 3.3.2 Robin boundary . . . . . . . . . . . . . . . 3.4 Space dependant diffusivity . . . . . . . . . . . . . 3.5 Discussion and future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 39 40 40 43 49 54 55 66 71 75 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 2.7 3 Testing . . . . . . . . . . . . . . . . Analysis . . . . . . . . . . . . . . . 2.6.1 Discrimination of tracks . . 2.6.2 Exogenous vs endogenous Discussion and future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendices A Additional derivations . . . . . . . . . . . . . . . . . . . . . A.1 From the Neumann splitting probability, equation 3.51 A.1.1 The Neumann Green’s function for the sphere . A.1.2 Derivation of χ1 . . . . . . . . . . . . . . . . . . A.2 From the Robin splitting probability, equation 3.70 . . A.2.1 Derivation of χ R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B Programs . . . . . . B.1 spotDetector3D B.2 mySpotDetector B.3 detecttr . . . . . B.4 testdetect . . . B.5 testtrack . . . . . . . . . . . . . . . . . 90 . 91 . 99 . 108 . 112 . 116 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 81 81 86 88 88 C Additional plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 v List of tables 2.1 2.2 2.3 2.4 2.5 Summary of the specifications of the experiments I analyzed data from . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Summary of the detection methods used . . . . . . . . . . . . Number of tracks which are effectively moving in each experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Maximum and minimum values of bootstrapped average velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Confidence intervals of 95% of vexo − vendo . . . . . . . . . . . 5 18 30 33 34 vi List of figures 1.1 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 2.13 2.14 2.15 2.16 2.17 Schematic diagram showing the major structural components of a lymph node . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Two photon microscopy requires two simultaneous events of high wavelength . . . . . . . . . . . . . . . . . . . . . . . . 6 Two photon microscopy presents less background noise . . . 6 A sample microscopy image from each data set . . . . . . . . 8 A sample image in original and after denoising in the two levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Sample trajectory computed by u-track . . . . . . . . . . . . 13 Extract from one of the movies to show how mySpotDetector works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 Resulting saved regions after processing the example shown in Figure 2.6 with mySpotDetector . . . . . . . . . . . . . . . 17 Two 3D views of the connected region from the example in Figure 2.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Sample simulated images . . . . . . . . . . . . . . . . . . . . 19 Error 1 as a function of the number of simulated cells . . . . 20 Error 2 as a function of the number of simulated cells . . . . 21 Histograms for an endogenous data set using mySpotDetectorbup as detection method with the lower level of denoising . 23 Histograms for an endogenous data set using u-track as method of detection with the lower level of denoising . . . . . . . . . 24 Histograms for an endogenous data set using u-track as method of detection with the higher level of denoising . . . . . . . . 25 Histograms for an exogenous using u-track as method of detection with the higher level of denoising . . . . . . . . . . . 26 Histogram of turning angle ignoring the z-coordinate for chosen tracks from all experiments with the higher level of denoising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3D representation of all tracks found for a data set . . . . . . 28 vii List of figures 2.18 Planar projections of Figure 2.17 . . . . . . . . . . . . . . . . 2.19 Turning angle histograms for odd and even time points separatedly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.20 Example of how the process often misses part of the trajectory 2.21 Histograms of step size corresponding to the endogenous experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.22 Histograms of step size corresponding to the exogenous experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.23 Distributions of bootstrapped average cell velocity . . . . . . 2.24 Cells change shape when they are on motion . . . . . . . . . 2.25 Microscopy images showing the contour of the LN . . . . . . 2.26 A trajectory found by eye which the computational process missed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Plots of the mean and variance of the FPT to the only target inner sphere with Neumann boundary condition . . . . . . . 3.2 Plot of the splitting probability to the only target inner sphere with Dirichlet boundary condition . . . . . . . . . . . . . . . 3.3 Plots of the mean and variance of the FPT to the only target inner sphere with Dirichlet boundary condition . . . . . . . 3.4 Plot of the splitting probability to the only target inner sphere with Robin boundary condition . . . . . . . . . . . . . . . . . 3.5 Plot of the splitting probability to the only target inner sphere with Robin boundary condition as a function of λ . . . . . . 3.6 Plots of the mean and variance of the FPT to the only target inner sphere with Robin boundary condition . . . . . . . . . 3.7 Plot of the mean and variance of the FPT to the only target inner sphere with Robin boundary condition as a function of λ 3.8 Model of the LN with 3000 APCs . . . . . . . . . . . . . . . . 3.9 Plots of the splitting probability to the one specific target along two radii when all targets have Dirichlet BC and the LN has Neumann BC . . . . . . . . . . . . . . . . . . . . . . . 3.10 Plots of the splitting probability to the one specific target along two circles when all targets have Dirichlet BC and the LN has Neumann BC . . . . . . . . . . . . . . . . . . . . . . . 3.11 Plots of the MFPT to the one specific target along two radii when all targets have Dirichlet BC and the LN has Neumann BC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 29 30 31 32 33 35 36 37 3.1 42 45 50 53 53 54 55 60 61 61 62 viii List of figures 3.12 Plots of the MFPT to the one specific target along two circles when all targets have Dirichlet BC and the LN has Neumann BC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Plots of the space-averaged MFPT to the one specific target as a function of the radius of the targets when all targets have Dirichlet BC and the LN has Neumann BC . . . . . . . . . . 3.14 Plots of the space-averaged MFPT to the one specific target as a function of ǫ and the number of targets when all targets have Dirichlet BC and the LN has Neumann BC . . . . . . . 3.15 Plots of the MFPT to the one specific target along two radii when all targets and the LN have Neumann BC . . . . . . . 3.16 Plots of the MFPT to the one specific target along two circles when all targets and the LN have Neumann BC . . . . . . . 3.17 Plots of the splitting probability to the one specific target along two radii when all targets have Dirichlet BC and the LN has Robin BC . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Plots of the splitting probability to the one specific target along two circles when all targets have Dirichlet BC and the LN has Robin BC . . . . . . . . . . . . . . . . . . . . . . . . . 3.19 Two-dimensional version of the FRN model . . . . . . . . . . 3.20 A possible function to describe how the diffusion coefficient changes on the FRN . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Comparison of the MFPT to the only target inner sphere with different BC on the LN . . . . . . . . . . . . . . . . . . . . . . C.1 Histograms for data set LN2 using mySpotDetector-bup as detection method with the lower level of denoising . . . . . C.2 Histograms for data set LN8 using mySpotDetector-bup as detection method with the lower level of denoising . . . . . C.3 Histograms for data set LN9 using mySpotDetector-bup as detection method with the lower level of denoising . . . . . C.4 Histograms for data set z1 using mySpotDetector-bup as detection method with the lower level of denoising . . . . . . . C.5 Histograms corresponding to data set LN10 assuming a pure Brownian motion model for the tracking step . . . . . . . . . C.6 Histograms corresponding to data set LN10 assuming a random motion plus movement with constant velocity model for the tracking step . . . . . . . . . . . . . . . . . . . . . . . . 62 63 63 66 67 70 70 71 73 75 119 120 121 122 123 124 ix List of figures C.7 Histograms corresponding to data set LN10 when the motion model assumed for the tracking step was random motion and movement along a straight line but with the possibility of immediate direction reversal . . . . . . . . . . . . . . C.8 Histograms corresponding to the analysis of the endogenous data set LN2 using u-track as method of detection and with the higher level of denoising . . . . . . . . . . . . . . . . . . . C.9 Histograms corresponding to the analysis of the endogenous data set LN8 using u-track as method of detection and with the higher level of denoising . . . . . . . . . . . . . . . . . . . C.10 Histograms corresponding to the analysis of the endogenous data set LN9 using u-track as method of detection and with the higher level of denoising . . . . . . . . . . . . . . . . . . . C.11 Histograms corresponding to the analysis of the endogenous data set Z1 using u-track as method of detection and with the higher level of denoising . . . . . . . . . . . . . . . . . . . . . 125 126 127 128 129 x List of abbreviations 2PM Two-photon microscopy APC Antigen presenting cell BC Boundary conditions COM Center of mass DC Dendritic cell FPT First passage time FRN Fibroblastic reticular network GFP Green flourescent protein HEV High endothelial venules LHS Left hand side LN Lymph node MAX Local intensity maxima MFPT Mean first passage time PLO Peripheral lymphoid organs RHS Right hand side SPT Single particle tracking TCR T cell receptor VFPT Variance of the first passage time xi Acknowledgements I would like to thank my advisor, Daniel Coombs, who took me as his student when I was desperate because none of the UBC faculty I wrote to had replied to me. I really enjoyed being his student and part of his group. A special thank you goes to Dodo Das who was always there to help. I feel his guidance was esential in the development of the first part of my thesis. Thanks to Isabell Graf who took the time to teach me the basics of homogenization theory, and to her advisor Peter Malte for listening to our ideas and giving his advice as a homogenization expert. I would also like to acknowledge the financial support and learning opportunities provided by the International Graduate Training Center (by the Pacific Institute for the Mathematical Sciences) throughout my masters, as well as the financial support granted to me by Consejo Nacional de Ciencia y Tecnologı́a from Mexico. Finally, I would like to thank the city of Vancouver and Canada for welcoming me, and all my friends and family, here and in Mexico, who helped me through the good and bad moments. xii Dedication To my beloved Mexico. xiii Chapter 1 Introduction 1.1 Biological description of the problem The immune system is protects the body from foreign agents and disease. In vertebrates the actions of the immune system can be classified in two main categories: innate immunity and adaptive immunity. The innate immune system acts first and provides a first line of defense. It is generic, meaning that it reacts in much the same way to similar pathogens. The adaptive immune system, on the other hand, identifies the pathogen and produces cells which will only attack that specific pathogen. The adaptive system keeps a copy of the cells created specifically to attack the pathogen even after it has been cleared out of the body. By doing this it provides long-lasting immunity to the host, so next time the same pathogen is found in the body the immune system will already have the specific tools to fight it. However, the initial adaptive response is slower and also relies on the innate system. T cells are a type of lymphocyte (white blood cell) which play a very important role in the adaptive immune response. T cells originate in the bone marrow but then migrate to the thymus to undergo maturation, which is why they are called “T” cells. Mature T cells then enter the bloodstream and migrate to the peripheral lymphoid organs (PLO), organized tissues where adaptive immune responses are initiated and where lymphocytes are maintained. T cells are continually recirculating through these tissues. They enter the PLOs by squeezing between high endothelial venules (HEV) and exit through efferent lymph vessels. Antigen1 is also carried from the site of infection to the PLOs. Specialized cells, referred to as Antigen Pressenting Cells (APC), carry antigen to the PLOs and display it to lymphocytes. [19] The lymph nodes (LN) are one type of PLOs. They are highly organized and are located at the points of convergence of vessels of the lymphatic sys1 Antigen - Term used to refer to any substance that can be recognized by the adaptive immune system. 1 1.2. Current approaches Figure 1.1: Schematic diagram showing the major structural components of a lymph node [24]. tem. One important structure inside the LN is the fibroblastic reticular network (FRN), a network formed by fibroblastic reticular cells and filaments between them. Figure 1.1 shows the architecture of the lymph node. A key step in the cascade of actions which lead to the adaptive response is T cell activation. T cells are equipped with receptors which allow them to recognize antigen displayed by APCs. Hence T cell receptors (TCR) show an enormous variability (from one cell to other). While in the lymph node, T cells successively scan APCs and read the information presented by them until finding their cognate antigen. If the reading contact with the cognate APC lasts long enough then the T cell is activated and it undergoes proliferation and differentiation, the first steps of adaptive immunity. The descendants of an activated T cell are ready to fight the pathogens (by activating B cells to produce antibody, or by directly killing infected cells). In particular, dendritic cells (DC) are an important class of APCs. They have long projections called dendrites, which gives them their name. 1.2 Current approaches While interacting with APCs, T cells change their behaviour in different ways depending on whether they are bound or not, or if they have already been activated. Therefore, analyzing the movement of T cells in the LN is important, for it can lead to better knowledge of the processes that take 2 1.2. Current approaches place inside the LN. Here, I present results from two different mathematical approaches to describing the behaviour of T cells in the LN. The first approach focused on describing T cell motion from data while the second one was about describing the time it takes for a T cell to find its corresponding APC under different hypotheses about T cell motion. The purpose of both approaches is to gain insight into the strategies of the immune system in making T cell activation an efficient process. In Chapter 2 I present the work done using the first approach. For this part I analyzed imaging data of labelled T cells inside the LN of a healthy mouse. Thanks to the development of advanced imaging techniques (twophoton microscopy) we were able to acquire time series of quasi-threedimensional images showing T cells inside a LN. From these I extracted cell paths and characterized them. This has been done in the past, however, we propose a new experimental frame which we think provides more accurate results. Here I present a comparison between the results obtained from the typical experiment and those from our new experimental frame. The second approach is presented on Chapter 3. Here we studied the first passage time of a model T cell to a fixed target, the cognate APC. We considered different scenarios by varying the motility conditions for the T cell as well as the characteristics of the environment where it moves, i.e. the LN. 3 Chapter 2 Analysis of two-photon imaging data from T cells moving within living tissue This chapter is the result of a collaboration with Dr. Pasquale Maffia’s team at the University of Glasgow. The data they provided us with was a time series of 3-dimensional images showing T cells inside the LN. Each 3D image corresponding to a time point is in reality a stack of several 2D images, spaced by a couple of micrometers in height, which are taken rapidly and sequentially, and so are regarded as forming one 3D image at a single time point. I worked with several data sets, which varied in the number of time points as wells as the size of the stack. Details of each data set’s specifications are given in Table 2.1. 2.1 Two-photon microscopy Two-photon excitation microscopy [6] is a relatively new fluorescence imaging technique. It has the advantage of allowing living tissue to be imaged at high depth. 2PM relies on the fact that two photons of approximatedly half the energy needed for single photon excitation, can equally excite a fluorophore if the two photons hit it simultaneously. Because that is a lower probability event, fewer fluorophores will become excited, and more excitation events are likely to occur near the focal point, where more photons are concentrated (Figure 2.1). This is in part the reason why 2PM presents less background noise (see Figure 2.2). Moreover, since lower energy means longer wavelength, the wavelength of the absorbed photons in 2PM is about two times the wavelength the absorbed photons in single photon microscopy. Because lower wavelength gets scattered a lot, this again means less noise in 2PM images, as well as less tissue damage by phototoxicity. [17] 4 LN2 LN7 LN8 LN9 LN10 Z1 Exogenous/ Endogenous endogenous endogenous endogenous endogenous exogenous exogenous Number of slices 17 12 12 30 22 25 Number of time points 120 150 60 20 30 38 Imaging rate (s/stack) 29.025 20.004 20.816 50.799 45.666 12.821 Pixels per side 512 512 512 512 512 560 µm per side 531.37 531.37 531.37 531.37 425.1 283.4 Height per slice (µm) 2 2 2 2 2 2.5 Table 2.1: Summary of the specifications of the experiments I analyzed data from. Only those specifications which influence the computational analysis of the data are included. 2.1. Two-photon microscopy Name 5 2.2. Experiments Figure 2.1: Two photon microscopy requires two simultaneous events of high wavelength [17]. Figure 2.2: Two photon microscopy presents less background noise [17]. 2.2 Experiments Previous studies of T cell motion are based on experiments where cells are extracted from the mouse, then exogenously tagged with a fluorophore protein, and finally transfected back into the mouse. In these experiments only some of the transfected cells then home to the LN. Hence images from this type of experiment show low densities of labelled cells. Moreover, we hypothesize that cells which have been outside the mice show altered dynamics, and also those cells rediscovered in the LN may represent an unusual subset of the basic T cell population. Our collaborators in the University of Glasgow have performed endogenous experiments on genetically modified mice whose cells express Green Fluorescent Protein (GFP). Because in these genetically modified mice all T cells express GFP, the density of cells shown in images obtained from endogenous experiments is higher. This is good in the sense that more data is available but it also makes the determination of trajectories more compli6 2.3. Methods for the inference of cell paths cated. 2.2.1 Experimental methods To image endogenous T cell populations, inguinal lymph nodes were excised from hCD2 GFP mice and placed in a perfusion chamber. Lymph nodes were maintained at 37°C by perfusion with heated CO2 independent media. Following staining, approximately 2 million T cells were transferred by intravenous tail vein injection into CD11c YFP transgenic mice. Multiphoton fluorescence was performed using a Zeiss 7MP multiphoton microscope (Carl Zeiss, Jena, Germany) after excitation by a Coherent Chameleon wavelength tuneable laser (Coherent, Glasgow, UK). Zeiss 10X/0.3 NA air and 20X/1.0 NA water immersion lenses (both Carl Zeiss, Germany) were used to focus light on the sample and to collect emitted light. Fluorescence emission was collected at <375nm (autofluorescence), 500-550nm (green) and 600-680nm (red). All samples were excited at 920nm. 2.3 Methods for the inference of cell paths 2.3.1 Single particle tracking Single particle tracking (SPT) is a widely used technique whereby the motion of biological particles can be observed. The particle of interest is tagged using any of a variety of experimental techniques. Then a series of images is taken by means of a microscope at high temporal resolution. From these one can attempt to extract particle trajectories, a highly non-trivial task that requires joining the dots (cell centers) between frames in a sensible way. In typical SPT experiments, only a small number of particles are tagged, which makes the tracking somewhat easier [18]. The endogenous experiments we worked with present three major difficulties as compared to typical SPT experiments: • First, in a strict sense, what I did was not SPT but rather single cell tracking. Cells are much bigger than what we generically call “a particle”. This is not really a disadvantage per se, in fact I would say it is an advantage. The problem, however, is that many existing SPT techniques cannot be directly applied to cell tracking. • Second, the high cell density is a big problem. Even by eye it is often difficult to tell which of two cells located close to each other on a 7 2.3. Methods for the inference of cell paths (a) LN2 (b) LN7 (c) LN8 (d) LN9 (e) LN10 (f) Z1 Figure 2.3: A sample image from each data set. Notice that the bottom row shows the exogenous experiments. All images correspond to t = 1 and z = 1. 8 2.3. Methods for the inference of cell paths frame moved to where on the next frame. Moreover, background noise increases when having so many bright spots on an image. • Third, SPT has only rarely been done in three dimensions. Having the third dimension is good because it provides more information about the phenomenon we are studying. However, we found that even once it is known which planes a specific cell lies on, determining its 3D center is a difficult task which has a great influence in the subsequent track analysis. Hence we can conclude that our problem is how to deal with so much information. 2.3.2 Steps in the tracking process I have distinguished three main stages in the process of extracting the particle trajectories. The stages were separated in this way because I tackled each stage separately, computationally speaking. Denoising This is the first stage and consists of filtering the images so that the background noise does not interfere with the subsequent stages. In other words, in this stage a noisy image is transformed into one with better defined cell boundaries. Detection The next step is to determine what makes a single cell in the image and find all single cells. To facilitate the last step once all cells in an image have been determined each of them is characterized by the 3D coordinates of its center. Tracking At this stage one has a list of 3D positions corresponding to each time point. The tracking step consists of connecting these positions into discrete trajectories in a sensible way. The following subsections describe each stage in more detail. A whole section (2.4) has been dedicated to the second stage, detection, because it was the step for which I developed new algorithms. One such algorithm actually uses the same method I used for the tracking step. This is why I left the detection section after the first and third steps’ subsections, which basically used methods that already existed. 9 2.3. Methods for the inference of cell paths 2.3.3 Denoising with wavelets As I said before, first the images need to be treated so as to get rid of noise and extract positions of objects that are truly there. Further, one has to dismiss all those objects which seem too small, or otherwise “flawed”, to be real T cells (or whatever cells we are tracking). For the denoising step I used an algorithm developed to extract bright spots from biological images [21]. The novelty about this algorithm is that it can detect bright spots at different levels of resolution, it is thus well suited to detect cells. The algorithm was implemented in MATLAB by Henry Jaqaman and François Aguet. The code includes both the denoising and the detection part. However, they both only work in 2D, which is why I split the code in two parts. The denoising part of the algorithm is based on image processing with wavelets theory. I thought about extending it to three dimensions but it is not clear how much that would help. The problem is that the resolution of the images is lower in the third dimension (height) than in each x-y plane. I thought if I used a direct extension to 3D of the denoising algorithm a lot of information might be discarded as noise. Hence, I decided to use the denoising part to process all 2D images individually. Then the 2D denoised images are put together as 3D images for the next step. At first I used the algorithm with the default parameter for the level of denoising and I noticed that after processing the images with this level of denoising there were many short tracks. So I decided to increase the level of denoising. The higher level of denoising uses more morphological operations to get rid of isolated pixels and other outliers. This is done by means of the MATLAB function bwmorph, the higher level additionally performing the processes erode, spur, clean, and thicken. By doing that I would get rid of dim spots in the central area, where most cells are located (see Figure 2.4). I tried both levels and most results here use the higher level of denoising, unless otherwise stated. However, there are advantages and disadvantages to both levels. The higher level sometimes misses spots that are clearly visible by eye in some frames breaking a path or splitting a cell in the z direction. 2.3.4 Track extraction U-track [13] is a well known SPT software. It is very flexible in that it allows the user to modify many parameters in order to adjust the algorithm to a 10 2.3. Methods for the inference of cell paths (a) Original (b) Lower level of denoising (c) Higher level of denoising Figure 2.4: An example image in original and after denoising in the two levels. This image corresponds to data set “LN2” with t = 1 and z = 8. specific problem. In short, what the algorithm does is the following: For each pair of consecutive time points, match positions in the first frame with positions in the following frame. This is treated as a linear assignment problem [2]: there is a cost for each match and the goal is to minimize the total cost. To reduce the size of the combinatorial problem, a position can only find a match in a circular area around it (the radius of which is a user-assigned parameter). Not having a match is highly penalized, but it may happen, as each particle can only find a single match. After the first step one is left with a bunch of generally short tracks. The second step is called “gap closing”, here the short tracks previously obtained are matched, this time not on consecutive time steps. At this stage u-track offers the possibility of allowing for merging and splitting. It often happens that two particles in an image come so close to each other that they look as a single cell. The concepts of merging and splitting are related to such events. Merging refers to two cells which had separated tracks before but at a given time point appear as a single position. Splitting, on the other hand, refers to the case when a single track splits into two different ones. If one allows for merging and splitting then two segments’ ends (starts) can be matched to the same segment’s start (end). The problem with using the merge & split capability is that at the end one has groups of tracks instead of a single cell track. For example, if two tracks merge and later split again into two separate tracks, there is no way of knowing which of the first two tracks corresponds with which of the two last tracks. For that reason in the results shown here I did not allow for merging and splitting. Next is a list of the u-track parameters I modified. 11 2.3. Methods for the inference of cell paths Gap closing time window This sets the maximum number of time steps that a track can be lost and still reappear. Because we have 3D images and the cells are large, bright objects, it is not often that a cell goes out of focus and then comes back. In fact, a cell can only get out of the imaging region if it is on the bottom or top planes, or near the edges of the other planes. Hence I set this parameter to 1. Minimum track segment length Tracks obtained from the first step which are shorter than this parameter (units of number of frames) will not be considered for the second step (the gap closing). For this parameter I used the default, which is 2. Search radius lower limit It refers to the search radius mentioned in the description of the first step of the tracking algorithm. Because many cells, specially in the periphery, stay static through out the video, I let this parameter be 0µm. Search radius upper limit This parameter was set to 10µm. Its value was determined by noting how big were the steps taken by some of the cells whose movement was more apparent by eye. I also considered the fact that the average diameter of a T cell is 7µm. Flag for linear motion U-track allows for three motion models: Brownian motion, Brownian motion plus movement with constant velocity, and random motion plus movement along a straight line but with the possibility of immediate direction reversal. I tested all three and the differences did not appear to be significant (see Figures C.5 to C.7 in Appendix C). The results shown in the main text were obtained using the third motion model. It is worth mentioning that during the detection step all calculations use units of pixels but before proceeding to the tracking step it is important to scale the coordinates appropiately so that the difference in x,y and z resolution does not introduce bias to the results. To that effect the results of detection were transformed to micrometers using the resolution specifications from each data set. 12 2.3. Methods for the inference of cell paths 19.5 19 18.5 18 17.5 17 158 16.5 156 16 154 15.5 120.5 120 119.5 119 118.5 152 118 117.5 117 116.5 Figure 2.5: Sample trajectory computed by u-track. This corresponds to LN2. The scale is micrometers. Notice that it spans about 6µm on each dimension. By comparing with the video one realizes that this cell did not really move. Its apparent movement is a consequence of a general drift by the whole lymph node. The zigzagging pattern is very likely a consequence of intensity variations in the cell throughout the video. 13 2.4. The detection step 2.4 The detection step 2.4.1 Spot detector The detection method from [21] (hereafter called spotDetector) is relatively simple so it was not hard to adapt it to work on three dimensions. Here is a description of the original spotDetector algorithm: 1. Find local intensity maxima (MAX). In the original 2D version of SpotDetector, a 9x9 window is used to define “local”. 2. Find connected regions. Two pixels are in the same region if they have one common vertex or side (hence each pixel has 8 neighbors). 3. For each connected region find its weighted centroid (COM) and the local maxima (computed in 1) that lies on the region. The program keeps a list of the COM of and highest local maximum in each region. 4. If no local maximum lies on the region then the index of the highest value on the region is recorded instead. 5. The program also keeps a list of all other (non highest) local maxima on the region, these are called secondary maxima. Secondary maxima whose distances to both the COM and the MAX do not excede a certain threshold (I used 5 pixels) are deleted from the list. 6. At the end, the COM and the remaining secondary maxima are returned as the centers of the detected cells comprised in the connected region or cluster. The program can easily be modified so the returned centers include the MAX instead of the COM. The only tricky part in making this algorithm work for the 3D data was defining local in the first step. The 2D version uses a MATLAB built-in function called ordfilt2 which replaces the value on each pixel with the maximum over the values of its neighboring pixels, where the neighborhood is defined by the notion of local (a 9x9 window). MATLAB does not have an analogous 3D function and although it does not seems like a complicated process to code, it turns out to be very computationally expensive. However, I found that someone had already coded a function ordfilt3 [5] which nevertheless has the restriction of only working on cubic windows to define locality. This was inconvenient because the x-y scale is very different from the z scale on my images. Nevertheless, I used the 3D version of spotDetector with ordfilt3 using a size 7 cubic window to define locality. Since scales 14 2.4. The detection step vary from data set to data set (see Table 2.1) this might be more suitable for some data sets than for other. One last detail from the 3D version is the definition of connected region. In 2D an 8-neighbor lattice was used. For 3D I chose 6-neighbors, i.e. two voxels are in the same region if they share a face. Like in the original spotDetector, the program offers the possibility of returning the COM or the MAX (plus the remaining secondary maxima). My 3D version of spotDetector, which has large parts from the original, is given in Appendix B.1. 2.4.2 My detection technique In the spotDetector algorithm, the detection of intensity local maxima and discrimination of them by distance is performed in order to distinguish the centers of several cells that may be very close together, forming a cluster. I proposed an alternative mechanism to separate cells using the 3D information. The idea behind my algorithm (hereafter referred to as mySpotDetector) is that clustered cells might appear as a single 2D connected region in some planes but not in all the planes spanned by all cells in the cluster. Hence one might take a 3D connected region and look for the planes in it where such region does not appear as a 2D connected region. The algorithm is explained next. 1. Find all 3D connected regions. I chose to do this defining connectedness in terms of a 6-neighbors matrix, but I also did trials using 18and 26-neighbors. 2. Loop over all these connected regions. 3. Look at the first (bottom most) plane that comprises the current region. 4. Find all 2D connected regions and save them. For these I used an 8neighbors definition (but also did trials with 4-neighbors definition). 5. Look at the following plane that comprises the current 3D region, and also find all 2D connected regions. 6. Look at the intersections between 2D regions in the current plane and saved regions. • If a region in the current plane intersects none of the saved regions, append it to the list of saved regions. 15 2.4. The detection step (a) z = 1 (b) z = 2 (c) z = 3 (d) z = 4 (e) z = 5 (f) z = 6 (g) z = 7 (h) z = 8 (i) z = 9 (j) z = 10 (k) z = 11 (l) z = 12 Figure 2.6: Extract from one of the movies to show how my detection technique works. All the planes contained in one three dimensional connected region are shown. Intensity variations have been omitted in this example. See next figure for the resulting saved regions. • If a region in the current plane intersects one or more saved regions, replace each of those saved regions with the corresponding intersection. 7. Repeat steps 5 and 6 until reaching the last plane that comprises the current 3D connected region. 8. A weighted centroid (depending on intensities) is then found for each saved region. These centers correspond to the cells that were so close to each other that appeared as the current single connected 3D region. 9. One could also start from the last (top most) plane comprising a given 3D connected region and sweep all planes down. This approach will be called “top down” whereas the one described before will be called 16 2.4. The detection step “bottom up”. It does not appear to make much difference which approach is used. Figures 2.6 to 2.8 show an example of how this detection method works. Figure 2.7: Resulting saved regions after processing the example shown in Figure 2.6 with mySpotDetector. According to step 8, the resulting detected cells would have centers determined as the intensity weighted centroids of the spots shown in this figure. Figure 2.8: Two 3D views of the connected region from the example in Figure 2.6. These were plotted with the Volume Viewer plug-in from ImageJ [20]. 2.4.3 Using u-track for detection Although the results from both previous methods did not seem too bad at first, when feeding them to u-track for the tracking step, the program would have difficulty processing the information and often would show errors in 17 2.4. The detection step computing the matrices to solve the linear assignment problem. Because of these errors we wanted another detection method. When consulting with the authors of u-track, it was suggested to us that we used u-track in 2D for the detection part. As it turned out, after doing this for the detection part the tracking step would then go smoothly. In this algorithm, u-track is used separately on the stacks corresponding to each time point. Here the z-position is the analogous of time in the normal tracking process. U-track parameters are set to allow only for very short motion and no gaps (because one cell cannot disappear and then reappear skipping one intermediate plane). The parameters used were: gap closing time window 1, minimum track segment length 5 (to reduce the number of tracks included in the gap closing step), search radius 1 (except for “z1” which has higher resolution and so the search radius was allowed to be 1 or 2), and the linear motion model was assumed to be random. The series of 2D positions input to u-track are obtained by using the original spotDetector on each plane individually. The output of using u-track as a detection method is a series of positions in the planes spanned by each cell. Coordinates of each cell are then determined in the following way: • For x and y, take the average of the x and y positions in each of the planes the spot spans. • The z-coordinate is also an average of the vertical positions of the planes where the spot appears. This average is weighted with the intensity of the cell on that plane normalized by the total intensity of the cell. Abbreviated name Program name sd3D-MAX sd3D-COM mySD-bup mySD-tdown u-track B.1. B.1. B.2. B.2. B.3. spotDetector3D spotDetector3D mySpotDetector mySpotDetector detecttr Section where it is described 2.4.1 2.4.1 2.4.2 2.4.2 2.4.3 Specific parameter choice = 0 choice = 1 bt = 0 bt = 1 Table 2.2: Summary of the detection methods used. See Appendix B for the codes to which the second column refers. The specific parameter is used in the codes to distinguish between two variations of the same algorithm. 18 2.5. Testing 2.5 Testing I mainly focused on testing the algorithm I created, that is mySpotDetector, although I did try a couple of tests on the existing methods (u-track and spotDetector). To test the detection methods I generated some 3-dimensional positions which were to be the centers of the cells. These centers are constrained to a box of size 254x254x64 and cells are assumed to be spherical of constant radius r. Then, the bottom plane of the box is “imaged”: an image corresponding to that plane shows in black those pixels which fall inside a sphere. The next three planes are not imaged, but the fourth is. 17 images are generated in this manner and then processed with the detection methods. See Figure 2.9 for sample simulated planes. The interested reader can find the code for this test on Appendix B.4. 50 50 50 100 100 100 150 150 150 200 200 250 200 250 50 100 150 200 250 250 50 (a) 10 100 150 200 250 50 (b) 100 50 50 50 100 100 150 150 150 200 200 200 250 250 100 150 (d) 1000 200 250 150 200 250 200 250 (c) 500 100 50 100 250 50 100 150 (e) 2500 200 250 50 100 150 (f) 5000 Figure 2.9: Sample simulated images with r = 6. Caption of each subfigure indicates number of simulated cells. It was difficult to come up with a way to measure error and this part still has to be improved. I computed two errors: Error 1 is the sum of the distance between each simulated cell and the closest center found by the detector, normalized by the number of cells simulated; Error 2 is the relative error between the estimated and true number of cells. Notice that Error 2 is not precisely the fraction of cells missed since many of the computed positions do not correspond to any of the real positions. 19 2.5. Testing Figure 2.10: Error 1 as a function of the number of simulated cells. The vertical axis units are micrometers and it is shown in log scale. It was not possible to test u-track when there are only a few or too many simulated cells (cases with 10,100, and 5000), since in such cases the algorithm returns no found cells. Figures 2.10 and 2.11 summarize the results obtained. The large relative error between the estimated and true number of cells (Error 2) is due to the high percentage of cell overlapping (see Figure 2.9). Taking that into account, one can conclude that the fact that Error 1 is so small, especially when there are too many simulated cells, is also due to cell overlapping. It implies that cluster centers are quite close to the true centers of the simulated cells. Notice that, at least in terms Error 2, the method that does best is mySpotDetector. The reason is that this method attempts to distinguish individual cells in a cluster by using only shape information, while the other methods rely on intensity variations (which the simulated images do not have). I tried to include the detection method which uses u-track in this testing procedure but I encountered several problems. When there are only a few cells or too many of them simulated (cases with 10,100, and 5000) the algorithm returns no found cells. For the cases where some cells were found, u-track returned several warning messages which indicate that it is not getting along well with the simulated data. So this method of detection 20 2.5. Testing Figure 2.11: Error 2 as a function of the number of simulated cells, shown as percentage. It was not possible to test u-track when there are only a few or too many simulated cells (cases with 10,100, and 5000), since in such cases the algorithm returns no found cells. does worse in the test than the other methods. The problem is that the 2D positions returned in the first step of the algorithm are not cell centers, but cluster centers. Clusters probably vary too much from one plane to the next one (in fact, mySpotDetector uses this variability to distinguish cells) and that is why u-track cannot build trajectories succesfully. I think this does not necessarily mean that the methods are bad, but maybe the simulated data is too far from its experimental counterpart and this is why they fail. Indeed, the failure comes from the fact that there are no local maxima since there is no variation in intensity at all. More tests are needed on simulated data that includes noise and different shapes for the cells. Even though u-track has been tested, results shown in [13] are only for 2D, so that it is not clear if the method was also tested in 3D. Therefore, I created a function (testtrack, see Appendix B.5) which generates a number of simulated Brownian paths with a chosen diffusion coefficient, saving their positions at certain time points. I wanted to build the simulation in such a way that positions were only recorded if the particle is inside the imaged box (i.e. a rectangular box of low height as before) but too much 21 2.6. Analysis information was lost and the tracker did a very bad job, returning only a handful of trajectories even when I was simulating 1000 tracks. This issue already raises a warning, although we had already suspected that the low resolution in the z direction was the biggest problem we had with the data. Nevertheless, I decided to go one step back and test u-track in a free diffusion situation, that is without the box constraints. Again I had trouble finding a good way of measuring the error. Just comparing the original number of tracks with those output by u-track does not seem to give a lot of information. In [13] the authors present their results from testing in terms of percentage of misses, but they do not explain how they assess whether a given trajectory was “missed” or not. In summary, testing needs improving, especially in terms of better measuring whether or not the results are acceptable. 2.6 Analysis From the results of the tracking process I analyzed some properties that can be rapidly computed and give an idea of the behaviour of the cells in terms of motion. These properties are Turning angle Angle formed by three positions of a particle in consecutive time frames, as measured by the dot product. Step size Distance travelled by a particle in one time step. Average step size Average step size over all steps taken by a given particle. Track diameter Maximum distance between any two positions of a given track. Track length (or life span) Number of time points a given track spans. Effectively computed as the difference between the last and the first time points where the particle is present. For each set of results I got I computed these properties for every trajectory and then made histograms showing this information. For example, Figures 2.12 and 2.14 both show results from data set LN7 but the first one corresponds to the analysis using the detection algorithm mySpotDetector while the second one corresponds to using u-track as method of detection. Notice from Figures 2.12(a), 2.13(a), and 2.14(a) that for the endogenous experiments there seem to be a preference for turning almost completely in 22 2.6. Analysis LN7sc − Turning angle 6000 2.5 Frequency Frequency x 10 3 5000 4000 3000 2 1.5 2000 1 1000 0.5 0 0 LN7sc − Step size 4 3.5 Pi/4 Pi/2 3Pi/4 0 0 Pi 20 Angle 40 (a) LN7sc − Average step size LN7sc − Diameter 3000 600 2500 Frequency 500 Frequency 80 (b) 700 400 300 2000 1500 1000 200 500 100 0 0 60 Velocity (um per minute) 10 20 30 40 0 0 20 um per minute 40 60 80 100 120 Distance (um) (c) (d) LN7sc − Track length 5000 Frequency 4000 3000 2000 1000 0 0 50 100 150 Time steps (e) Figure 2.12: Histograms for data set LN7 using mySpotDetector-bup as detection method with the lower level of denoising. 23 2.6. Analysis LN7trd − Step size LN7trd − Turning angle 1200 100 1000 Frequency Frequency 80 60 40 800 600 400 20 0 200 Pi/4 Pi/2 0 0 3Pi/4 5 (a) LN7trd − Diameter 20 LN7trd − Track length 1000 200 800 Frequency Frequency 15 (b) 250 150 100 50 0 0 10 Distance (um) Angle 600 400 200 5 10 Distance (um) (c) 15 20 0 0 20 40 60 80 100 Time steps (d) Figure 2.13: These histograms correspond to the analysis of the endogenous data set LN7, using u-track as method of detection and with the lower level of denoising. 24 2.6. Analysis the opposite direction (compare also with Figure 2.15 which shows results for one of the endogenous experiment). The peak at an angle of π would not be expected for random motion. In fact, no preferred angle would be expected for random motion. LN7trd − Turning angle LN7trd − Step size 800 80 700 70 600 Frequency Frequency 60 50 40 30 500 400 300 20 200 10 100 0 Pi/4 Pi/2 0 0 3Pi/4 5 Angle 10 15 20 Distance (um) (a) (b) LN7trd − Diameter LN7trd − Track length 200 700 600 500 Frequency Frequency 150 100 400 300 200 50 100 0 0 5 10 15 20 Distance (um) (c) 0 0 20 40 60 80 100 Time steps (d) Figure 2.14: These histograms correspond to the analysis of the endogenous data set LN7, using u-track as method of detection and with the higher level of denoising. 2.6.1 Discrimination of tracks The mentioned bias in turning angle could be an important feature but we first needed to rule out the possibility that it was a consequence of the methods used and not of the true intrinsic cell dynamics. On Figure 2.16 25 2.6. Analysis LN10trd − Turning angle LN10trd − Step size 25 300 250 Frequency Frequency 20 15 10 200 150 100 5 50 0 Pi/4 Pi/2 0 0 3Pi/4 2 4 Angle 6 8 10 12 Distance (um) (a) (b) LN10trd − Diameter LN10trd − Average step size 100 80 70 80 Frequency Frequency 60 50 40 30 20 60 40 20 10 0 0 2 4 6 8 0 0 10 10 20 30 40 Distance (um) Distance (um) (c) (d) LN10trd − Track length 160 140 Frequency 120 100 80 60 40 20 0 0 5 10 15 20 25 30 Time steps (e) Figure 2.15: These histograms correspond to the analysis of exogenous data set LN10, using u-track as method of detection and with the higher level of denoising. 26 2.6. Analysis LN2trd − Turning angle ignoring z−coordinate LN7trd − Turning angle ignoring z−coordinate 50 50 40 Frequency Frequency 40 30 20 10 0 0 30 20 10 Pi/4 Pi/2 3Pi/4 0 0 Pi 3Pi/4 Pi LN8trd − Turning angle ignoring z−coordinate LN9trd − Turning angle ignoring z−coordinate 2 8 Frequency 1.5 6 4 0 0 1 0.5 2 Pi/4 Pi/2 3Pi/4 0 0 Pi Pi/4 Pi/2 3Pi/4 Pi Angle Angle LN10trd − Turning angle ignoring z−coordinate Z1trd − Turning angle ignoring z−coordinate 25 6 5 Frequency 20 Frequency Pi/2 Angle 10 Frequency Pi/4 Angle 15 10 4 3 2 5 0 0 1 Pi/4 Pi/2 Angle 3Pi/4 Pi 0 0 Pi/4 Pi/2 3Pi/4 Pi Angle Figure 2.16: Histogram of turning angle ignoring the z-coordinate for all experiments with the higher level of denoising. Only tracks with average step > 3µm and with at least 5 frames are included in these histograms. Compare the endogenous versus the exogenous experiments (bottom row). 27 2.6. Analysis we ignored the z-coordinate information. We did so because given the way this coordinate is assigned, there is a lot of oscillation in it, specially when the cell is only taking short steps (see Figures 2.17 and 2.18). 35 30 25 20 15 10 5 0 600 600 400 500 400 200 300 200 0 100 Figure 2.17: 3D representation of all tracks found for LN2 (higher level of denoising, u-track as method of detection). Scale is in micrometers. 600 500 35 35 30 30 25 25 20 20 15 15 10 10 400 300 200 100 0 100 5 150 200 250 300 350 400 (a) x - y 450 500 550 0 100 5 150 200 250 300 350 400 (b) x - z 450 500 550 0 0 100 200 300 400 500 600 (c) y - z Figure 2.18: Planar projections of the tracks extracted from LN2. Scale is in micrometers. Another approach was to separate the list of positions from a trajectory into two lists: odd and even time points. Then compute the histograms for each list separately. If the cell, or one of its coordinates, was oscillating in such a way that it was going back and forth from one “side” to another, then these histograms would show the true turning angle distribution. This 28 2.6. Analysis LN8trd − Turning angle skipping even positions LN8trd − Turning angle skipping odd positions 35 25 30 20 Frequency Frequency 25 20 15 15 10 10 5 5 0 Pi/4 Pi/2 Angle 3Pi/4 0 Pi/4 Pi/2 3Pi/4 Angle Figure 2.19: Turning angle histograms for odd and even time points separatedly. This results correspond to LN8 processed with the higher level of denoising and u-track as tracking method. The z-coordinate, unlike stated in Section 2.4.3, was assigned as the simple average of the planes spanned by a cell (i.e. no weighting). was not the case however, as it can be seen from Figure 2.19. Because we observed many very short tracks (see insets (d) and (e) from Figure 2.15, and bear in mind that the diameter of a T cell is about 7µm) we decided to restrict our analysis to the longer tracks. We also thought that the bias in turning angle could be being introduced by the short tracks as well. Moreover, I noticed for a few examples that the average distance travelled during a time step is less than 1µm for tracks that correspond to stationary cells. For that reason I computed the number of tracks whose average step is greater than 3µm/min, which is roughly equivalent to more than 1µm per frame. Results are shown on Table 2.3 as well as on Figure 2.16. By comparing the tracking results with the videos I concluded that most of the short tracks are part of a bigger trajectory that the algorithm failed to connect (see Figure 2.20). Hence the fact there are so many short tracks does not mean, as one could imagine, that cells move so much that they disappear from the imaged region in just a few frames. 29 2.6. Analysis Name LN2 LN2 LN7 LN7 LN8 LN8 LN9 LN9 LN10 LN10 Z1 Z1 Level of denoising higher normal higher normal higher normal higher normal higher normal higher normal Tracks w/ big avg steps 540 753 326 530 161 248 23 43 98 111 76 158 Long tracks 100 140 56 90 30 34 3 9 40 41 14 23 Table 2.3: Number of tracks which are effectively moving, according with the criterion of having average step greater than 3µm/min. The third column shows how many of the moving tracks are at least 5 steps long. Shorter tracks do not contribute significant information. (a) (b) Figure 2.20: (a) Example of how the process often misses part of the trajectory. The blue line shows what the computer found to be the full track, the magenta line is the rest of the track, which I found by eye. For each time point in the trajectory I chose the z plane where the cell looked bigger in area and then add all those images up to obtain the trace shown in the picture. The extract comes from experiment LN2. (b) 3D view of the part of the track which was detected by u-track. 30 2.6. Analysis LN2trd − Step size LN7trd − Step size 6000 1000 800 4000 Frequency Frequency 5000 3000 600 400 2000 200 1000 0 0 10 20 30 40 50 0 0 60 10 Velocity (um per minute) LN8trd − Step size 600 600 500 Frequency 500 Frequency 30 40 50 60 50 60 LN9trd − Step size 700 400 300 400 300 200 200 100 100 0 0 20 Velocity (um per minute) 10 20 30 40 50 60 Velocity (um per minute) 0 0 10 20 30 40 Velocity (um per minute) Figure 2.21: Histograms of step size corresponding to the endogenous experiments. 2.6.2 Exogenous vs endogenous We thought that an interesting feature to compare between the exogenous and endogenous experiments was the cell velocity. Two of the properties described on Section 2.6 are related to cell velocity, namely step size and average step size. Step size refers to the distance traveled by a given cell in a given single time step, whereas average time step is the mean of the step sizes of all steps from a given track. Both quantities can be scaled to measure distance traveled per minute (which has actually been done in all the plots displayed here). I used the step size information to test our hypothesis that cells present altered behaviour in exogenous experiments. We asked ourselves if, given the experimental results that we have, was there enough information to reject the hypothesis that the average velocity is the same for both types 31 2.6. Analysis LN10trd − Step size Z1trd − Step size 350 20 300 15 Frequency Frequency 250 200 150 100 10 5 50 0 0 10 20 30 40 50 60 0 0 Velocity (um per minute) 50 100 150 200 Velocity (um per minute) Figure 2.22: Histograms of step size corresponding to the exogenous experiments. The obvious outliers from experiment Z1 were deleted before bootstrapping. (endogenous and exogenous) experiments. I found that we can reject that hypothesis with a confidence level of at least 95%. To take advantage of all the data that we have, we used step size instead of average step size. In this way I have 17,261 data points from the endogenous experiments and 1,349 data points from the exogenous experiments. To obtain the distribution of mean velocity from this data, I took 1000 bootstrap samples from each type of experiment and computed the average of those samples. The bootstrap sampling was done in MATLAB. The hypotheses to test are H0 : vexo − vendo = 0 Ha : vexo − vendo > 0 where vexo refers to the average velocity of a cell from an exogenous experiment while vendo refers to the average velocity of a cell from an endogenous experiment Because the cells from one of the exogenous experiments, namely LN10, exhibit very rapid movement as compared to the others, I bootstrapped the data from the two exogenous experiments separately as well as all together. A high amount of movement is a natural consequence of the experimental conditions which is in general reduced by other methods but can be worse in certain experiments, like it is the case for this particular movie. Tables 2.4 and 2.5 and Figure 2.23 show the results obtained. From Table 2.5 we can see that cells from the exogenous experiments move significantly 32 2.6. Analysis faster than those from endogenous experiments. In fact, from Table 2.4 we see that the distributions of average velocity do not overlap, so technically our conclusion has a significance level of 0%. Although we would need to analyze more experiments before drawing such a strong conclusion, this result is highly suggestive of an important difference in labelled T cell motion between exogenous and endogenous experiments. Similar results were obtained when using the normal level of denoising. Variable vendo vexo vZ1 v LN10 Minimum 2.0902 3.7452 11.0362 2.5323 Maximum 2.2662 5.0438 16.1753 2.9707 Table 2.4: Maximum and minimum values of bootstrapped average velocities. Units are µm/min. vendo is the bootstrapped average of 1000 samples collected from all four endogenous experiments while vexo is the corresponding bootstrapped average of 1000 samples from the two exogenous experiments. Because the cells from LN10 exhibit very rapid movement as compared to the others, the data from the two exogenous experiments was bootstrapped separately as well as all together. 0.07 0.06 Probability 0.05 0.04 0.03 0.02 0.01 0 2 4 6 8 10 12 14 16 Velocity (um per minute) Figure 2.23: Distributions of bootstrapped average cell velocity: vendo in red, vexo in black, vZ1 in magenta, and v LN10 in blue. 33 2.7. Discussion and future work Variable vexo − vendo vZ1 − vendo v LN10 − vendo Confidence interval 1.6772 - 2.3442 9.3529 - 12.9622 0.4314 - 0.7080 Table 2.5: Confidence intervals of 95% of vexo − vendo . Units are µm/min. 2.7 Discussion and future work Although I did not obtain as many results from this project as I would have liked, I gained a lot of insight into the computational problem it poses. I think the reason why tracks are apparently easy to identify by eye is because our eyes and brain incorporate a lot of information, including shape and history. Unless a cell has too many neighbors, one can often determine unequivocally what is the position of a cell in the next frame by indirectly using information about its shape and its past trajectory. Hence, shape plays an important role in this problem. It is true that all cells have more or less the same shape, and so it is not possible to perform the tracking by identifying the cells like one would with people. However, when cells move they change their cytoskeleton following more or less the same patterns which might be helpful in the tracking process. On the other hand, static cells look more round (see Figure 2.24). I think one ought to consider, if not the shape, at least the whole area (or volume) of a cell instead of just attempting to form tracks out of a list of centers. From one time point to the next one the area occupied by a cell in the first frame often overlaps the area it occupies in the following frame. I am certain that by considering the whole volume of the cell and not just its center, the problems brought by the oscillations in coordinates, such as the strong peak in turning angle at π, would be overcome. There is another piece of information that I think should be incorporated into u-track: the fact that cells do not simply disappear from one time step to the next. Because we are imaging in 3D they do not go out of the imaging region unless they are in the bottom/top planes or near the edges of the planar images. I predict that this would prevent the algorithm from breaking trajectories into several pieces, at least partially. One other problem I did not manage to solve is the general drift that the whole lymph node undergoes. It is evident from the movies, when one looks at the static cells, that there is a general drift going on. Such drift is not obvious when the movie is played but only when one looks at the first 34 2.7. Discussion and future work (a) t = 6 (b) t = 7 (c) t = 8 (d) t = 9 Figure 2.24: These images correspond to frames 6 to 9, plane z = 1, of experiment LN2. Notice how the shapes of the marked cells change from one frame to another. The cell marked with a blue dot is static while the one marked with a red dot is not. and last frames. It is because of this drift that the diameters and average steps from the tracks corresponding to static cells are not so small. For all the T cell images we have, we also have images where the boundary of the LN is tagged (Figure 2.25). We thought that from this images we were going to be able to estimate the size and direction of the drift. I tried to do so by computing the intensity weighted centroid of the LN’s boundary for each time point, and then compute the vector of movement of this centroid from one frame to the next one. This vector would then be substracted from all cell positions in all future frames. I did this with one set but in doing so I noticed that some cells actually seem to be drifting in a very different direction than what the vector of weighted centroid displacement showed. Finding a better way to measure the drift is on the list of future work. These and other factors are the reason why we think the tracking results need a lot of improvement yet. Lack of adequate visualization tools is one more problem that, if overcome, could help improve the results. For example, at one point I engaged in the task of generating a considerable amount of trajectories by eye. I ended up not using these results simply because they were written on a piece of paper. I am convinced that it would have been different if the software we used to play the movies, which is only in 2D, had allowed me to store as part of the movie the trajectories I was detecting. What I was planning to do with these trajectories was train the parameters of the tracking, and maybe also the detection, algorithms. Visualization of the movies with the superimposed tracks in 3D is also important to validate the results and detect problems in the methods. I do not think I have the right training 35 2.7. Discussion and future work Figure 2.25: Microscopy images showing the contour of the LN. These images correspond to the first and last time points, plane z = 1, from data set LN7. 36 2.7. Discussion and future work Figure 2.26: This figure shows a trajectory found by eye which the computational process missed. For each time point in the trajectory I chose the z plane where the cell looked bigger in area and then add all those images up to obtain the trace shown in the picture. The extract comes from experiment LN2. to build my own software to solve this problem. However, I spend some time looking for existing solutions and even talked to some experts on visualization of biological images. Unfortunately, I did not find a satisfactory solution. Besides improving the detection and tracking methods one natural extension to this work is to analyze more experiments, especially of the exogenous type, in order to confirm the results from Section 2.6.2. Further extensions of this work include implementing some statistical method to infer the trajectory of a cell between time steps or once it has left the imaging region. Also, evaluating the relationship between movement patterns and zones in the LN, the influence of zones being not only suggested by theory but also by observing the videos. 37 Chapter 3 First passage time of a T cell to an APC In this chapter I present some calculations of the mean time it takes for a naı̈ve T cell to reach a specific APC, among many non-specific ones, inside the LN. First passage time calculations have been performed in the past for other problems related to T cell activation [27]. The T cell is taken to be a point in three dimensional space. The APCs, which in reality are bigger than T cells, are model as spheres, first with fixed finite radii and then with asymptotically small volume. While T cells are assumed to be mobile, APCs are static. The lymph node is in most cases a sphere, for the sake of simplicity in explicit calculations, but most results hold for any sort of convex shape. The idea of analyzing this problem came from two papers [11, 12] which studied a similar problem using computational tools. The authors simulated the movement of T cells in the LN, including other cells in the model. On the first paper the goal was to find how population-level parameters relate to single-cell properties in the process of pathogen killing by T cells. It was found that the time it takes for a T cell to kill an infected cell has a bigger influence on the killing efficency than the time needed for a T cell to find the target cell. The second paper studied the influence of the FRN in the probability of two cells finding each other. The results in this chapter are theoretical, that is, not based on nor compared with data. Hence, a predetermined model has to be assumed for cell motion. For the most part, I assumed that T cells perform isotropic Brownian motion. In Section 3.4 I considered the case of space dependant diffusion. Another factor involved in the calculations is the influence of the lymph node’s boundary. Because T cells enter the lymph node mainly through HEVs, located along veins that cross the LN, they could show up inside it at virtually any point. They however would always exit through the same spot, which could be modelled as a hole on the LN’s boundary. Nevertheless, cells do not exit when they first “bump into” this hole, but they con38 3.1. First passage time concepts tinue their search until they receive a signal to exit [24]. Hence I considered various different scenarios for the boundary conditions on the LN. Since T cells are assumed to perform Brownian motion, then their concentration u ( x, t) at time t on the point x in space, satisfies the diffusion equation ∂ ∂u ∂u = D ( x) ∂t ∂x ∂x The different boundary conditions considered apply to this equation and are the following Dirichlet u ( x, t) = 0 for x ∈ ∂Ω, which means the particles get absorbed by the boundary (or exit) as they reach it. Neumann ∂n u ( x, t) = 0 for x ∈ ∂Ω, where ∂n represents the derivative in the outward normal direction, means that there is no outward flux through the boundary. In other words, when a particle hits the boundary it is reflected back inside Ω. Robin The Robin boundary condition is a combination of the Dirichlet and Neumann conditions, namely ∂n u ( x, t) = −λu ( x, t) for x ∈ ∂Ω. With this boundary condition there is some flux through the boundary. 3.1 First passage time concepts An important concept in the theory of stochastic processes is that of the first passage time of a certain event. In the case of Brownian motion we usually think of the first time the particle reaches certain target. The probability moments of the first passage time are of interest in many areas of study [22]. In the next sections I will constantly refer to the following two concepts: Splitting Probability When the absorbing (Dirichlet) boundary consists of N B , P is the probability of exiting through B . disjoint sets B = ∪i=1 i i i Moments of the first passage time T n is the n-th moment of the first passage time to the target(s) of interest, that is the APC. When referring to the mean first passage time (MFPT), that is the first moment, we drop the subscript and say simply T. 39 3.2. Basic problem 3.2 Basic problem On a first approach to the problem we consider Ω, the LN, to be the interior of a sphere of radius R+ , and Ω1 , an APC, to be a concentric sphere of radius R− < R+ . For simplicity we center the spheres at the origin. For this problem, we considered two different conditions for the boundary of Ω. The boundary of the target, ∂Ω1 , is always absorbing, as we regard the process as terminated once the particle has reached Ω1 . The T cell is a point with initial position x ∈ ΩΩ1 . 3.2.1 Neumann boundary Although in reality T cells are not strictly trapped inside the LN, as I already said, it is also not true that they exit as soon as they hit the boundary of the LN. Therefore, one might be interested in asking, if the T cell was indeed unable to leave the LN, how long will it take it on average, to find its cognate APC. Mean first passage time It is well known [22, 23] that the MFPT T satisfies the Poisson equation 1 , x ∈ ΩΩ1 ; D T (x) = 0, x ∈ ∂Ω1 ; ∂n T (x) = 0, x ∈ ∂Ω, ∇2 T (x) = − (3.1) where D is the diffusion constant. Because of radial symmetry T only depends on r so that the partial differential equation 3.1 reduces to the ordinary differential equation 2 ∂T 1 ∂2 T (r) + (r) = − , r ∈ ( R− , R+ ), 2 ∂r r ∂r D T (R− ) = 0, ∂n T (R+ ) = 0. (3.2) The solution to the homogeneous version of the Poisson equation, that is the Laplace equation, is c1 + c2 , (3.3) r 40 3.2. Basic problem 2 r and we find a particular solution for 3.2 to be − 6D . Hence the general solution to 3.2 is r2 c1 + c2 − . (3.4) T (r ) = r 6D Now the boundary conditions T ( R− ) = 0 and ∂n T ( R+ ) = 0 give the following equations for c1 and c2 c1 + c2 − R− c1 T ′ ( R+ ) = − 2 − R+ R2− = 0, 6D R+ = 0, 3D from where c1 = − c2 = R3+ , 3D R2 R3 + 2R3+ R2− c R3+ − 1 = −+ = − . 6D R− 6D 3DR− 6DR− The final expression for T is T (r ) = R3− + 2R3+ r2 R3 . − + − 6DR− 3Dr 6D (3.5) A plot of T is shown on Figure 3.1. Variance of the mean first passage time Since the equation for all the moments of the FPT are also well known, it is possible to also obtain an expresion for the Variance of the First Passage Time (VFPT) using that Var = T2 − T2 . (3.6) As in [23], the equation for the second moment is 2 T (x) , x ∈ Ω; D T2 (x) = 0, x ∈ ΩΩ1 ; ∂n T2 (x) = 0, x ∈ ∂Ω. ∇2 T2 (x) = − (3.7) I found the following to be a particular solution for the previous equation: R3− + 2R3+ r2 R3+ r r4 + − . 60D2 3D2 18D2 R− 41 3.2. Basic problem Figure 3.1: Plot of the mean (left) and variance (right) of the FPT to the target inner sphere Ω1 as a function of the starting point, with reflecting boundary conditions. R+ = 1, R− = 0.01, and D = 3. Hence the solution to 3.7 is R3− + 2R3+ r2 r4 R3+ r c1 + c2 + + − . T 2 (r ) = r 60D2 3D2 18D2 R− After substituting the boundary conditions and solving for c1 and c2 we find a final expression for the second moment of the FPT 18R5+ R− − 5R3+ R3− + 2R3+ T2 (r) = 45D2 R− r r 3R− r3 + 60R3+ R− − 10r R3− + 2R3+ (3.8) + 180D2 R− 40R6+ + 7R6− − 72R− R5+ − 20R3− R3+ , + 180D2 R2− and so we also have an expression for the VFPT, which we denote by V r 20R3+ − r3 10R6+ + R6− − 20R3− R3+ − 36R− R5+ + V (r) = 90D2 90D2 R2− (3.9) 2R5+ R6+ + − . 5D2 r 9D2 r2 A plot of V is shown on Figure 3.1. 42 3.2. Basic problem 3.2.2 Dirichlet boundary Next we study what would happen if T cells would actually exit the LN as soon as they hit its boundary. If we only changed the Neumann boundary conditions to Dirichlet BC on equation 3.1 we would be computing the mean time it takes for the particle to be absorbed, which now can happen on any of the boundaries, i.e. ∂Ω or ∂Ω1 . However, that is not what we want. We want the MFPT to the APC, that is the mean time it takes for the particle to hit ∂Ω1 before hitting the LN’s boundary ∂Ω. Hence, as in [22], we first need to know the probability of absorption at Ω1 , that is the splitting probability P1 . Because we will only compute the splitting probability to Ω1 we drop the subscript 1 and use P instead of P1 throughout this chapter. Splitting probability and mean first passage time According to [22] the splitting probability and the MFPT to Ω1 satisfy the following equations: ∇2 P (x) = 0, x ∈ ΩΩ1 , P (x) = 0, x ∈ ∂Ω, P (x) = 1, x ∈ ∂Ω1 ; D ∇2 [P T] (x) = −P (x) , x ∈ ΩΩ1 , [P T] (x) = 0, x ∈ ∂Ω ∪ ∂Ω1 . (3.10) (3.11) Again by spherical symmetry the previous equations are reduced to 2 P ′′ (r) + P ′ (r) = 0, r P (R− ) = 1, P (R+ ) = 0; (3.12) 2 ′ ′ D P T + 2P T + P T + P T + PT = −P , r P (R− ) T (R− ) = 0, P (R+ ) T (R+ ) = 0. (3.13) ′′ ′ ′ ′′ 43 3.2. Basic problem We already now the form of the solution to 3.12 (equation 3.3) so we just need to check boundary conditions and get the constants. c1 + c2 = 1 R− c P ( R + ) = 1 + c2 = 0 R+ c1 =⇒ c2 = − R+ =⇒ R+ c1 − R− c1 = R− R+ R − R+ =⇒ c1 = R+ − R − R− =⇒ c2 = − R+ − R − P ( R− ) = Hence P (r ) = R− R − R+ − . r ( R+ − R − ) R+ − R − (3.14) Equation 3.13 then reduces to T ′′ + and we have that 1 2c2 T′ = − c1 + c2 r D (3.15) 2c2 2 =− , c1 + c2 r R+ − r so the solution is T (r ) = − r (r − 2R+ ) k1 − + k2 . 6D r − R+ Now we need to find k1 and k2 from the boundary conditions. We have P (r) T (r) = c2 R+ (r − 2R+ ) c2 r (r − 2R+ ) − 6D 6D k 1 c2 R + k 1 c2 k c R + − − 2 2 + + k 2 c2 r (r − R+ ) r − R+ r where we have written c1 in terms of c2 using that c1 = − R+ c2 . Notice that even though P ( R+ ) = 0, P ( R+ ) T ( R+ ) could be nonzero since T has a r −1R+ term. We have 0 = P ( R+ ) T ( R+ ) = k 1 c2 ( R + − R + ) k 1 c2 =− R+ ( R+ − R+ ) R+ 44 3.2. Basic problem Figure 3.2: Typical shape of the probability of reaching Ω1 , P , as given by equation 3.14. The parameters used are R− = 0.01, R+ = 1, and D = 3. from where we must have k1 = 0. On the other hand since P ( R− ) = 1 then to satisfy P ( R− ) T ( R− ) = 0 we need 0 = T ( R− ) = − and thus k2 = R− ( R− − 2R+ ) + k2 , 6D R− ( R− − 2R+ ) . 6D The final expresion for T is r (r − 2R+ ) R− ( R− − 2R+ ) + . (3.16) 6D 6D Plots of P and T for the case with absorbing boundary are shown on Figures 3.2 and 3.3. T (r) = − Variance of the first passage time In [22], the author does not present equations for the higher order moments of the FPT. However, it is easy to deduce the equation for the second moment following the same arguments used to obtain the equation for the MFPT [see 22, Section 1.6.3]. We obtain P ( x) by summing the probabilities for all paths that start at x and reach R− , i.e. (3.17) P ( x) = ∑ Pp ( x) . p∈ paths 45 3.2. Basic problem Then, by their definitions, T and T2 can be written as T ( x) = ∑ p∈ paths Pp ( x) t p ( x) , ∑ p∈ paths Pp ( x) (3.18) T2 ( x) = ∑ p∈ paths Pp ( x) t2p ( x) , ∑ p∈ paths Pp ( x) (3.19) where t p ( x) represents the time it takes for the particle to go from x to R− following path p. Now the sum overall the considered paths can be decomposed into the outcome after one step and the sum over all path remainders from the intermediate point x′ . In this discretizing framework we assume that on a time step δt the particle can only take steps of size δx in one of six directions (±e1 , ±e2 , ±e3 ), and each direction is chosen with equal probability. Thus x′ can only be equal to x ± δxei with i ∈ {1, 2, 3} . Therefore # " 3 3 1 1 P (x) = ∑ ∑ Pp (x + δxei ) + ∑ Pp (x − δxei ) 6 6 p i=1 i=1 = 1 3 P p (x + δxei ) + P p (x − δxei ) . ∑ 6 i=1 Upon substituting P ( x ± δxei ) by a three term Taylor expansion one obtains equation 3.10. If we repeat what we just did to equation 3.17, to equation 3.19 we obtain 46 3.2. Basic problem P (x) T2 (x) = Pp (x) t2p (x) ∑ p∈ paths 3 = ∑∑ p i=1 2 1 Pp (x + δxei ) t p (x + δxei ) + δt 6 δt 3 = δt P (x) + ∑ 3 i=1 2 " 2 1 + Pp (x − δxei ) t p (x − δxei ) + δt 6 ∑ Pp (x + δxei ) t p (x + δxei ) p + ∑ Pp (x − δxei ) t p (x − δxei ) p + 1 3 6∑ i=1 " # ∑ Pp (x + δxei ) t2p (x + δxei ) p + ∑ Pp (x − p δxei ) t2p (x − δxei ) # = δt2 P (x) δt 3 + ∑ [P (x + δxei ) T (x + δxei ) + P (x − δxei ) T (x − δxei )] 3 i=1 + 1 3 [P (x + δxei ) T2 (x + δxei ) + P (x − δxei ) T2 (x − δxei )] 6∑ i=1 where the last expression is obtained using equations 3.18 and 3.19. Next we replace [P T n ] ( x ± δxei ) by its three term Taylor expansion about x. Then 2 δt 3 2 ∂ PT 2 P (x) T2 (x) = δt P (x) + ∑ 2P (x) T (x) + δx (x) 3 i=1 ∂x2i δx2 ∂2 P T2 1 3 + ∑ P (x) T2 (x) + (x) 3 i=1 2 ∂x2i = δt2 P (x) + 2δtP (x) T (x) + + P (x) T2 (x) + δtδx2 2 ∇ [P T] (x) 3 δx2 2 ∇ [P T2 ] (x) 2 =⇒ −2δtP (x) T (x) δtδx2 2 δx2 2 = δt2 P (x) + ∇ [P T] (x) + ∇ [P T2 ] (x) . 3 2 (3.20) 47 3.2. Basic problem Now in this discretized diffusion process a relation between the time and space steps is enforced, namely that D= δx2 . 2δt Thus equation 3.20 turns into δx2 2D 2 P (x) + δx2 2 δx2 δx4 2 ∇ [P T] (x) + ∇ [P T2 ] (x) = − P (x) T (x) . 6D 2 D Finally, we divide the previous equation by δx2 and then take the limit as δx goes to 0 to obtain ∇2 [P T2 ] ( x) = − 2 P ( x) T ( x) . D (3.21) Boundary conditions are analogous to those in 3.13. Now that we have derived the equation we proceed to solving it. As in equation 3.15 the previous equation reduces to T2′′ + 1 2 T′ = (r (r − 2R+ ) − R− ( R− − 2R+ )) r − R+ 2 3D2 (3.22) which has solution given by T2 (r) = 3D2 b1 + 5R3+ R− (R− − 2R+ ) + 2R5+ 45D2 (r − R+ ) (R − R+ )2 (r − R+ )2 (r − R+ )4 − − + + b2 , 9D2 20D2 where b1 and b2 are constants to be determined upon sustitution of boundary conditions. We start by checking the condition imposed on the outer boundary (r = R+ ), again we must be careful here since T2 has a term r −1R− . Cancelling terms where applicable we are left with the following equation at r = R+ 0 = [P T2 ] (r) R+ 3D2 b1 + 5R3+ R− (R− − 2R+ ) + 2R5+ − r 3D2 b1 + 5R3+ R− (R− − 2R+ ) + 2R5+ = 45D2 r (r − R+ ) (r − R+ ) =− 3D2 b1 + 5R3+ R− (R− − 2R+ ) + 2R5+ . 2 45D r (r − R+ ) 48 3.2. Basic problem Hence 5R3+ R− ( R− − 2R+ ) + 2R5+ . 3D2 Now, since P ( R− ) = 1 at the other boundary we have b1 = − 0 = T2 (R− ) (R− − R+ )4 (R− − R+ )4 + + b2 =− 9D2 20D2 and thus 11 ( R− − R+ )4 . 180D2 Finally, we arrive to an expression for T2 11 ( R− − R+ )2 − 9 (r − R+ )2 ( R− − R+ )2 − (r − R+ )2 T 2 (r ) = . 180D2 (3.23) We just need to use this together with formula 3.6 and equation 3.16 to obtain the VFPT for the absorbing case: b2 = V (r) = 11 (R− − R+ )2 − 9 (r − R+ )2 = 180D2 2 (R− − R+ )2 − (r − R+ )2 (R− − R+ )2 − (r − R+ )2 − 36D2 2 (R− − R+ ) − (r − R+ )2 3 (R− − R+ )2 − 2 (r − R+ )2 90D2 (3.24) . Figure 3.3 shows the VFPT for the problem with Dirichlet boundary conditions. 3.2.3 Robin boundary Finally, we examine the case of a Robin boundary condition for ∂Ω. It is a good condition to consider since the T cell is not really unable to leave the LN, as in the reflecting case, but it does not necessarily exits the first time it hits the boundary, as in the absorbing case. Therefore, it makes sense to ask what would happen if, on reaching the boundary of the LN, the T cell would exit only with certain probability smaller than 1. With this assumption and following the same discretizing argument that was used 49 3.2. Basic problem Figure 3.3: Plot of the mean (left) and variance (right) of the FPT to the target inner sphere Ω1 as a function of the starting point, with absorbing boundary conditions (equations 3.16 and 3.24). The parameters used are R− = 0.01, R+ = 1, and D = 3. on the previous section, we derived a boundary condition for the splitting probability. As it turns out, this is in agreement with assuming a Robin condition for the diffusion equation for the particle concentration u (as in the introduction fo this chapter) [8]. The same equations 3.10, 3.11, and 3.21 are still valid, only the boundary conditions will change. To derive the proper boundary conditions we go back to one dimension, the result is easily extendable to our 3-dimensional case. We consider a particle moving on the interval [R− , R+ ]. Derivation of boundary conditions Suppose that if the particle reaches R+ then it either takes a step of size δx back into the interval with probability 1 − λ (δx), or the process ends with probability λ ( δx). In other words, the probability of getting to R+ − δx from R+ in a time step of size δt, p ( R+ − δx, δt| R+ ), is 1 − λ ( δx). Assume also that this probability function λ satisfies λ (δx) = λ0 δx, where λ0 is a constant. Hence λ (δx) → 0 when δx → 0. Then P (R+ ) = (1 − λ) P (R+ − δx) ≈ (1 − λ) P (R+ ) − δxP ′ (R+ ) , (3.25) 50 3.2. Basic problem where we have used a Taylor approximation to obtain the second line. Therefore (1 − λ) δxP ′ ( R+ ) = −λP ( R+ ) . Now we divide by δx and obtain (1 − λ0 δx) P ′ ( R+ ) = −λ0 P ( R+ ) . Taking the limit as δx goes to 0 we get the following condition for the splitting probability on the right hand side boundary: P ′ ( R + ) = − λ0 P ( R + ) . (3.26) Following analogous steps, an equation can be derived for the MFPT to R− , T ( x). From equation 3.18 and by our boundary assumption we have [P T] (R+ ) = ∑ (1 − λ) Pp (R+ − δx) t p (R+ − δx) + δt p = δt (1 − λ) P (R+ − δx) + (1 − λ) P (R+ − δx) T (R + − δx) ′ = δtP (R+ ) + (1 − λ) [P T] (R+ ) − δx [P T] (R+ ) , where we have subsituted equation 3.25 and again used a Taylor approximation of P T about R+ to obtain the last line. This yields (1 − λ) δx [P T ]′ ( R+ ) = δtP ( R+ ) − λ [P T ] ( R+ ) Recall that we have established a relation between δt and δx, namely δt = δx2 /2D. Substituting that and the definition of λ (δx) on the previous equation and then dividing by δx we get δx P ( R+ ) − λ0 [P T ] ( R+ ) . 2D Again taking the limit as δx goes to 0 gives the boundary condition we were looking for: (3.27) [P T ]′ ( R+ ) = −λ0 [P T ] ( R+ ) . (1 − λ0 δx) [P T ]′ ( R+ ) = The same condition is satisfied by the second moment, as can be seen following the analogous argument, i.e. 2 [P T2 ] (R+ ) = ∑ (1 − λ) Pp (R+ − δx) t p (R+ − δx) + δt p from where = δt2 P (R+ ) + 2δt (1 − λ) [P T2 ] (R+ − δx) + (1 − λ) [P T2 ] (R+ ) − δx [P T2 ] ′ (R+ ) [P T2 ]′ ( R+ ) = −λ0 [P T2 ] ( R+ ) . (3.28) 51 3.2. Basic problem Splitting probability and mean first passage time Now we are ready to solve the sytem given by equations 3.10, 3.11, and 3.21 changing the boundary conditions on Ω1 to conditions 3.26, 3.27, and 3.28 respectively, which corresponds to the FPT problem with a partially absorbing boundary. We drop the subscript on parameter λ0 in the boundary conditions. We first look at the equation for P . We already know the general solution to 3.10, namely c c P (r ) = 1 + 1 − 1 , r R− Then condition 3.26 gives c1 c1 c1 +λ 1− = 0, − 2 +λ R+ R+ R− (3.29) which then implies −c1 R− + λc1 R+ R− + λR2+ R− − λc1 R2+ = 0 and so c1 = c2 = λR2+ R− , λR2+ − λR+ R− + R− (3.30) R− (1 − λR+ ) . − λR+ R− + R− (3.31) λR2+ Hence the solution for P is R− λR2+ − λR+ r + r . P (r ) = r (λR2+ − λR+ R− + R− ) (3.32) See figures 3.4 and 3.5 for a visual representation of this equation. We can now work the equation for T. Again we know the equation simplifies by spherical symmetry to equation 3.15, which substituting the values of c1 and c2 we just computed yields T ′′ + 1 2 (1 − λR+ ) T′ = − . 2 λR+ − λR+ r + r D The general solution of the previous equation is T (r ) = k 2 − λR2+ − λR+ r + r 6D (λR+ − 1) 2 2 − k1 . (3.33) 6D (λR+ − 1) (λR2+ − λR+ r + r) 52 3.2. Basic problem Figure 3.4: Solution for the splitting probability with a Robin BC on ∂Ω. Parameters as before, plus λ = 0.5. Figure 3.5: Splitting probability P as a function of the “leakyness” of the boundary, λ, for a fixed value of the starting point, r0 = 0.7. All other parameters remain the same. The spectrum of λ has been split in two plots since P has a singularity at λ = 1 for the chosen value of R+ . 53 3.3. Multiple targets with asymptotically null volume Figure 3.6: Solution for the MFPT with a Robin BC on ∂Ω. Parameters as before, plus λ = 0.5. After substituting both boundary conditions, T ( R− ) = 0 and [P T ] ′ ( R+ ) = −λ [P T ] ( R+ ), we obtain an expression for the MFPT: " 2 2 λR2+ − λR+ R− + R− − λR2+ − λR+ r + r 1 T (r) = 3D (λR+ − 1) 2 (λR+ − 1) # R3+ (r − R− ) . − λR2+ − λR+ R− + R− λR2+ − λR+ r + r (3.34) See figures 3.6 and 3.7 for a visual representation of this equation. 3.3 Multiple targets with asymptotically null volume The real problem we are concerned with is how long it takes for a T cell to find its antigen specific APC, when it is surrounded by other non-specific APCs. The best model would be one where, upon reaching the non-specific APCs, the T cell in question would remain attached to it only for a short time and then it would continue its search in the LN. We start in this section with a model where the T cell stops its search when it hits any of the APCs. We consider N spherical targets Ωǫi with centers xi ∈ Ω with radii ǫai , ai > 0, in such a way that these spheres are completely contained inside 54 3.3. Multiple targets with asymptotically null volume Figure 3.7: MFPT T as a function of the “leakyness” of the boundary, λ, for a fixed value of the starting point, r0 = 0.7. All other parameters remain the same. The spectrum of λ has been split in two plots since T has a singularity at λ = 1 for the chosen value of R+ . the LN (a sphere Ω of radius a centered at the origin). There are two reasons for making the radii asymptotically small. First, for simplicity, since having randomly located spheres with positive volumes would break the spherical symmetry of Ω that we have been taking a lot of advantage from. By keeping the symmetry we can use a matched asymptotics rather than a simulations approach to solve the problem. The second reason, is that, although APCs tend to be larger than T cells, it makes more sense to make APCs have asymptotically null rather than positive volume, given that the T cell is merely a point. We again examine different boundary conditions for the LN. 3.3.1 Neumann boundary This problem has already been studied in the past [3, 4, 15]. In this section I will outline the process of obtaining a two term expansion of the splitting probability P1 , that is the probability of reaching Ωǫ1 before reaching any of the other APCs Ωǫ j for j = 2, . . . , N, by using matched asymptotics. The N Ω . I followed notes from starting point x is in Ω\Ω p where Ω p ≡ ∪i=1 ǫi Professor Michael Ward [25, 26]. For simplicity of notation, we again drop the subscript 1 on the splitting probability and use only P ( x). Later we 55 3.3. Multiple targets with asymptotically null volume will use numeric subscripts again for the asymptotic expansions, then P1 will be used again, but it should not be confused with what we are now simply denoting by P . Again P satisfies Laplace’s equation ∇2 P (x) = 0, ∂n P (x) = 0, P (x) = 1, P (x) = 0, x x x x ∈ Ω\Ω p ; ∈ ∂Ω; ∈ ∂Ω1 ; ∈ ∪N j=2 ∂Ωǫi . (3.35) We first look at the outer region, away from the targets. As usual, we expand P as P = P 0 + ǫ P 1 + ǫ2 P 2 + . . . (3.36) Then ∇2 Pk = 0, x ∈ Ω\ { x1 , . . . , x N } ∂n Pk = 0, x ∈ ∂Ω (3.37) with certain singularity conditions as x → x j for j = 1, . . . , N, determined upon matching to the inner solution. Then we look at the inner regions, we analyze separately the region near each target. Near the j-th APC, we expand the inner solution w (y) ≡ P x j + ǫy , with y ≡ ǫ−1 x − x j , w = w0 + ǫw1 + . . . (3.38) Now ∇2y w = ǫ2 ∇2x P and so ∇2y w0 = 0, y ∈ / Ω j ; w0 = δj1 , y ∈ ∂Ω j ∇2y w1 = 0, y ∈ / Ωj ; w1 = 0, y ∈ ∂Ω j (3.39) (3.40) Here Ω j = ǫ−1 Ωǫ j . The far-field boundary conditions for w0 and w1 are determined by the matching condition as x → x j between the inner and outer expansions 3.36 and 3.38, written as P0 + ǫP1 + ǫ2 P2 + . . . ∼ w0 + ǫw1 + . . . (3.41) The first matching condition is that w0 ∼ P0 as |y| → ∞. As before, the solution to 3.39 is c1 + c2 , w0 = |y| 56 3.3. Multiple targets with asymptotically null volume so that as |y| → ∞, w0 = c2 and thus P0 = c2 is a constant. The condition w0 = δj1 , y ∈ ∂Ω j , means w0 a j = δj1 , and so c1 = a j δj1 − P0 . Hence a j δj1 − P0 w0 ( y ) ≈ P0 + |y| (3.42) as |y| → ∞. From this and the matching condition 3.41 we conclude that as x → xj , a j δj1 − P0 P1 ∼ . (3.43) |x − xj | To be able to satisfy this condition we slightly perturb equation 3.37: ∇2 P1 ( x) = Sδ ( x) , where S is a constant and δ is the Dirac delta function. Transforming the previous equation to spherical coordinates we have 1 ∂ 2 ∂ P1 r = Sδ (r) , r2 ∂r ∂r which upon integration holds ∂ 2 ∂ P1 r = Sr2 δ (r) . ∂r ∂r (3.44) Now recall that by definition Z R3 δ ( x) dx = 1. We can transform this equation to spherical coordinates to obtain 1= = = = Z 2π Z π Z ∞ 0 0 Z ∞ 0 Z ∞ 0 Z ∞ 0 = 4π 0 δ (r) r2 δ (r) r2 δ (r) r Z ∞ 0 2 δ (r) r2 sin θdrdθdϕ Z 2π Z π 0 Z 2π 0 Z 2π 0 0 sin θdθdϕdr [− cos θ]0π dϕdr (3.45) 2dϕdr δ (r) r2 dr. 57 3.3. Multiple targets with asymptotically null volume Hence, using 3.45 on equation 3.44 we get ∂ P1 ∂r ∂ P1 ∂r S + c1′ 4π c1′ S =⇒ = + 4πr2 r2 S + c1′′ + c2′ . =⇒ P1 (r) = − 4πr r2 = Now we can find out what S has to be to agree with the matching condition 3.43, namely S = −4π δj1 − P0 a j . Therefore the equation for P1 writes N ∇2 P1 = −4π ∑ δj1 − P0 a j δ x − x j , x ∈ Ω; (3.46) j=1 ∂n P1 = 0, x ∈ ∂Ω. Next we will use the divergence theorem to state a solvability condition for P1 and thus determine P0 . Indeed by the divergence theorem Z ZZ Ω ∇2 P1 dV = ZZ ∂Ω ∇P1 ·~ndS = ZZ ∂Ω ∂n P1 dS and the RHS is equal to zero by the boundary condition on 3.46. On the other hand ZZZ N ∇2 P1 dV = −4π ∑ δj1 − P0 a j Ω j=1 and so ( 1 − P 0 ) a1 − P 0 ∑ N j=2 a j = 0 a =⇒ P0 = N 1 . ∑ j=1 a j We write P0 = a1 + · · · + a N a1 , ā ≡ , N ā N (3.47) (3.48) to save space. 58 3.3. Multiple targets with asymptotically null volume Now let us look at the Green’s function for the system 1 − δ (x − ξ) , x ∈ Ω; |Ω| ∂n G = 0, x ∈ ∂Ω; ∇2 G = Z Ω (3.49) G (x; ξ) dx = 0. In terms of this function the solution to 3.46 can be written as N P1 ( x) = 4π ∑ δj1 − P0 a j G x; x j + χ1 (3.50) j=1 where χ1 is a constant. By integrating the previous equation we obtain an expression for χ1 : Z N Ω P1 dx = 4π ∑ δj1 − P0 a j j=1 = χ1 | Ω | Z Ω G x; x j + Z Ω χ1 by equation 3.47. Then 1 χ1 = |Ω| Z Ω P1 dx We have finally arrived to a two term expansion of the splitting probability: " # 1 N a1 + 4πǫa1 G ( x; x1 ) − a j G x; x j + ǫχ1 + O ǫ2 . (3.51) P ( x) = ∑ N ā N ā j=1 The derivation of G and χ1 can be found in references [15, 25, 26] and I have also rewritten it in Appendix A for the sake of completeness. See Figures 3.9 and 3.10 for some example plots of P . Mean first passage time We now look at the MFPT. The solution to 1 , x ∈ Ω\Ω p , D ∂n T G (x) = 0, x ∈ ∂Ω, T G (x) = 0, x ∈ ∂Ω p , ∇2 T G (x) = − (3.52) 59 3.3. Multiple targets with asymptotically null volume 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 1 1 0.5 0.5 0 0 −0.5 −0.5 −1 −1 Figure 3.8: Sphere of radius 1 representing the lymph node, with 3000 randomly placed spherical targets (which represent dendritic cells). The one target filled in green is Ω1 , which represents the cognate APC the T cell is looking for. All targets are placed in such a way that there is at least 5ǫ between each pair and also separating them from the LN boundary. Since a spherical LN in mice is about 2mm in diameter [14], the scale here is in millimeters. I used a j = 0.35 for all j, since for ǫ = 0.01 the radius ǫa j = 3.5 × 10−3 mm would correspond to the average size of a dendritic cell. However, such a big value for ǫ made it imposible to compute positions for the 3000 targets without overlapping so I used ǫ = 0.001 instead. See the following figures for more details. 60 3.3. Multiple targets with asymptotically null volume −3 1.2 −3 x 10 2.5 1 x 10 2 0.8 P P 1.5 0.6 1 0.4 0.5 0.2 0 0 0.2 0.4 0.6 0.8 0 0 1 0.2 r 0.4 0.6 0.8 1 r Figure 3.9: Splitting probability from equations 3.51, A.11, and A.15, with ǫ = 0.001. P is plotted as a function of r along the radius shown in red (left) or black (right) in figure 3.8 (the angles which define these radii were chosen randomly). −4 5.5 −4 x 10 8 5 x 10 7 4.5 6 4 P P 5 3.5 4 3 3 2.5 2 2 1.5 0 1.5708 3.1416 Azimuthal angle 4.7124 6.2832 1 0 1.5708 3.1416 4.7124 6.2832 Azimuthal angle Figure 3.10: Splitting probability from equations 3.51, A.11, and A.15, with ǫ = 0.001. P is plotted as a function of the azimuthal angle along the circle shown in red (left) or black (right) in figure 3.8 (the radius and the polar angle which define these circles were chosen randomly). 61 3.3. Multiple targets with asymptotically null volume corresponds to the mean first passage time to hitting any of the traps. (The subscript G stands for “general”.) This is readily calculated by using the same matched asymptotic approach as before. Thus T G is given in the limit ǫ → 0 of small trap radius by " # N 4πǫ T a3 2 TG ≈ 1 − 4πǫ ∑ a j G x; x j + (3.53) a Ga + O ǫ 3N āDǫ N ā j=1 [4]. See Figures 3.11 and 3.12 for examples. 0.1 0.115 0.098 0.11 0.096 0.105 T T 0.094 0.1 0.092 0.095 0.09 0.09 0.088 0.086 0 0.2 0.4 0.6 0.8 0.085 0 1 0.2 r 0.4 0.6 0.8 1 r Figure 3.11: MFPT from equations 3.51 and A.11, with ǫ = 0.001 and D = 3. T G is plotted as a function of r along the radius shown in red (left) or black (right) in figure 3.8 (the angles which define these radii were chosen randomly). 0.1035 0.106 0.103 0.105 0.1025 T T 0.104 0.102 0.103 0.1015 0.102 0.101 0.101 0.1005 0.1 0 1.5708 3.1416 4.7124 Azimuthal angle 6.2832 0.1 0 1.5708 3.1416 4.7124 6.2832 Azimuthal angle Figure 3.12: MFPT from equations 3.51 and A.11, with ǫ = 0.001 and D = 3. T G is plotted as a function of the azimuthal angle along the circle shown in red (left) or black (right) in figure 3.8 (the radius and the polar angle which define these circles were chosen randomly). 62 3.3. Multiple targets with asymptotically null volume Also, the average R MFPT based on a uniform distribution of starting − 1 points, T̄ G = |Ω| Ω T G dx, satisfies T̄ G ≈ 4πǫ T a3 a G a + O ǫ2 . 1+ 3N āDǫ N ā (3.54) T̄ G for parameters like those from Figure 3.11 is around 0.1075. Increasing the radii of the targets to 3.5 gives a T̄ G of about 0.0125, while reducing ǫ to 0.0001 makes T̄ G = 1.0885. See also Figures 3.13 and 3.14. 0.35 4 3.5 0.02 0.3 0.018 3 0.25 0.016 2.5 0.2 0.014 0.15 0.012 2 1.5 1 0.008 0.05 0.5 0 0.01 0.1 0.02 0.04 0.06 0.08 0.1 0 0.006 0.2 0.4 0.6 0.8 1 0.004 2 4 6 8 10 Figure 3.13: Plot of T̄ G versus a j . N = 3000, a = 1, ǫ = 0.001 and D = 3. Due to the big variation in scale I separated the range of a j in three plots. 1.8 1 1.6 1.4 0.8 1.2 1 0.6 0.8 0.6 0.4 0.4 0.2 0.2 2 4 ε 6 8 10 0 500 1000 1500 2000 2500 3000 3500 4000 −4 x 10 N Figure 3.14: Plot of T̄ G versus ǫ and N. Same parameters as previous figure. To compute the MFPT to the specific cognate APC, Ω1 , when the T cell is unable to leave the LN we would instead need to solve 1 P (x) , x ∈ Ω\Ω p ; D ∂n T (x) = 0, x ∈ ∂Ω; [P T] (x) = 0, x ∈ ∂ ∪ N j=2 Ωǫ j ; ∇2 [P T] (x) = − (3.55) T (x) = 0, x ∈ ∂Ωǫ1 . 63 3.3. Multiple targets with asymptotically null volume Unfortunately, we have not as yet been able to solve this. Instead, we study the previously mentioned case where the T cell keeps its search after having contact with other non-cognate APCs. This is effectively done by imposing a Neumann boundary condition on all other targets (j = 2, . . . , N). Under this model, a contact between the T cell and a non-cognate APC is infinitesimaly short. The differential equation to solve is 1 , x ∈ Ω\Ω p ; D S ∂n T (x) = 0, x ∈ ∂Ω ∪ N j=2 Ωǫ j ; ∇2 T (x) = − (3.56) T (x) = 0, x ∈ ∂Ωǫ1 . In the outer region we expand T = ǫ−1 T0 + T1 + ǫT2 + . . .. For k = 0, 2 the T k satisfy ∇2 T k (x) = 0, x ∈ Ω\ { x1 , . . . , x N }, ∂T k (x) = 0, x ∈ ∂Ω, (3.57) while T1 satisfies 1 , x ∈ Ω \ { x1 , . . . , x N } ; D ∂n T1 (x) = 0, x ∈ ∂Ω. ∇2 T1 (x) = − (3.58) Singularity conditions as x → x j will be determined upon matching to the inner solution. In the inner region near the j-th trap, we let like we did for P before, − 1 − 1 v ( y) ≡ T x j + ǫy with y ≡ ǫ x − x j . We also expand v = ǫ v0 + v1 + . . .. The equations for the vk are ∇2 vk y = 0, y ∈ / Ωj , ∂n vk y = 0, y ∈ ∂Ω j , j 6= 1, (3.59) vk y = 0, y ∈ ∂Ω1 , and the matching condition reads ǫ−1 T0 + T1 + ǫT2 + . . . ∼ ǫ−1 v0 + v1 + . . . By equations 3.57 and 3.59, T0 is a constant and v0 is too in the inner region near the j-th trap for j 6= 1. The matching condition implies that v0 = T0 in the j-th inner region j 6= 1 while v0 ( y ) = T 0 − T 0 a1 |y| 64 3.3. Multiple targets with asymptotically null volume in the inner region near the first trap. By checking the matching condition again, we conclude that T1 satisfies equation 3.58 with singular behaviour T1 ∼ − a1 T 0 | x − x1 | as x → x1 and T1 ∼ 0 as x → x j . Therefore, in terms of the Dirac distribution, T1 satisfies ∇2 T1 (x) = 4πa1 T0 δ (x − x1 ) − ∂n T1 (x) = 0, x ∈ ∂Ω. 1 , x ∈ Ω; D We use the divergence theorem on the previous equation to derive a solvability condition for T1 , which yields T0 = |Ω| . 4πa1 D (3.60) Now in terms of the Green’s function, defined by equation 3.49, and a constant χ T we find that T1 = −4πa1 T0 G ( x; x1 ) + χ T . Next we expand the previous expression for T1 as x → x j using the local behaviour of G found on Appendix A.1.1. We get a1 T 0 − 4πa1 T0 R1,1 + χ T , as x → x1 ; | x − x1 | T1 ∼ −4πa1 T0 Gj,1 + χ T , as x → x j . T1 ∼ − Upon substituting into the matching condition we find that v1 ∼ −4πa1 T0 G j,1 + χ T , as |y| → ∞, where we are using the matrix notation from Appendix A.1.1. From equation 3.59 we note that v1 is also aconstant in the inner region near the j-th trap for j 6= 1, whereas v1 (y) = c1 1 − trap, with c1 a constant. Hence a1 |y| v1 (y) = χ T − 4πa1 T0 G j,1 in the inner region near the first δj,1 a1 1− |y| . 65 3.3. Multiple targets with asymptotically null volume 317.4 317.4 317.35 317.3 317.3 317.2 317.25 317.1 T 317 T 317.2 317.15 316.9 317.1 316.8 317.05 316.7 317 316.6 0 0.2 0.4 0.6 0.8 1 316.5 0 0.2 r 0.4 0.6 0.8 1 r Figure 3.15: MFPT from equation 3.61 with ǫ = 0.001 and D = 3. T is plotted as a function of r along the radius shown in red (left) or black (right) in figure 3.8 (the angles which define these radii were chosen randomly). Once again we go back to the matching condition and obtain that T2 should also be approximatedly 0 as x → x j , j 6= 1, and T2 ∼ a1 (4πa1 T0 R1,1 − χ T ) | x − x1 | as x → x1 . Thus the equation for T2 is ∇2 T2 (x) = −4π 4πa1 T0 R1,1 − χ T a1 δ (x − x1 ) , x ∈ Ω; ∂n T2 (x) = 0, x ∈ ∂Ω. Finally, the divergence theorem yields a solvability condition for T2 which gives χ T = 4πa1 T0 R1,1 . We have thus obtained a two-term expression for the MFPT: T ( x) = a3 4πa3 + ( R1,1 − G ( x; x1 )) . 3ǫa1 D 3D (3.61) See Figures 3.15 and 3.16 for examples. 3.3.2 Robin boundary In this section we want to study the asymptotic solution to 3.35 when the boundary condition on ∂Ω, ∂n P = 0, is replaced by ∂n P = λP . Since the Robin Green’s function for the Laplace equation is unknown, we instead consider the case when λ = ǫλ0 for some constant λ0 and let ǫ go to 0. Notice that the solution in the inner region would be the same from the previous section, as only the outer boundary condition has changed. 66 3.3. Multiple targets with asymptotically null volume 317.34 317.3 317.32 317.28 317.3 317.26 317.24 T T 317.28 317.26 317.22 317.2 317.18 317.24 317.16 317.22 317.2 0 317.14 1.5708 3.1416 4.7124 317.12 0 6.2832 1.5708 Azimuthal angle 3.1416 4.7124 6.2832 Azimuthal angle Figure 3.16: MFPT from equation 3.61 with ǫ = 0.001 and D = 3. T is plotted as a function of the azimuthal angle along the circle shown in red (left) or black (right) in figure 3.8 (the radius and the polar angle which define these circles were chosen randomly). Therefore, we need only solve the outer problem and be careful that the match is not lost. Like before, in the outer region we expand P = P0 + ǫP1 + ǫ2 P2 + . . .. Notice that the boundary condition for P0 is again ∂n P0 = 0 so that P0 is a constant. Now, near the j-th trap we again let w (y) ≡ P x j + ǫy and so equations 3.39 and 3.40 hold. From 3.43 we find that P1 must satisfy N ∇2 P1 = −4π ∑ δj1 − P0 a j δ x − x j , x ∈ Ω, (3.62) j=1 ∂n P1 = −λ0 P0 , x ∈ ∂Ω. The solvability condition for P1 obtained by the divergence theorem gives N −4π ∑ δj1 − P0 a j = j =1 ZZ ∂Ω λ0 P0 (3.63) = −λ0 P0 |∂Ω|, where |∂Ω| is the surface area of the sphere Ω. Thus P0 = 4πa1 . 4πN ā + λ0 |∂Ω| (3.64) 67 3.3. Multiple targets with asymptotically null volume Next we decompose P1 = −λ0 P0 P1P + P1H + χ R where χ R is a constant and P1P and P1H obey the following conditions |∂Ω| , x ∈ Ω, |Ω| = 1, x ∈ ∂Ω, ∇2 P1P = Z ∇2 P1H = Z ∂n P1P Ω (3.65) P1P (x) dx = 0; N λ0 P0 |∂Ω| − 4π ∑ δj1 − P0 a j δ x − x j , x ∈ Ω, |Ω| j=1 ∂n P1H = 0, x ∈ ∂Ω, Ω (3.66) P1H (x) dx = 0. In terms of the Green’s function of the previous section (equation 3.49), a solution to 3.66 is readily found to be N P1H = 4π ∑ δj1 − P0 a j G x; x j , (3.67) j=1 by identity 3.63. As for P1P , the divergence theorem yields Z Ω Z G (ξ; x) ∇2 P1P (ξ) − P1P (ξ) ∇2 G (ξ; x) dξ = ∂Ω G (ξ; x) dξ. The LHS is equal to λ0 P0 |∂Ω| |Ω| Z Ω G (ξ; x) dξ − Z Ω P1P (ξ) Hence P1P ( x) = Z ∂Ω 1 − δ (ξ − x) dξ = P1P (x) . |Ω| G (ξ; x) dξ. (3.68) Indeed, by reciprocity of G the third condition for P1P is satisfied too since Z Z Z G (ξ; x) dx dξ P1P (x) dx = ∂Ω Ω Ω Z Z G (x; ξ) dx dξ = ∂Ω Ω = 0. 68 3.3. Multiple targets with asymptotically null volume For the specific case we are dealing with, namely Ω is a sphere of radius a, we can also find a explicit solution to 3.65. First notice that |∂Ω| 4πa2 3 = 4 = . 3 |Ω| a 3 πa Thus the general solution for P1P is P1P (r) = c1 r2 + c2 + r 2a where c1 and c2 are constants. Then the boundary condition yields c c r 1 ′ |r=a = − 12 + 1 1 = P1P ( a) = − 2 + r a a so that c1 = 0. On the other hand, the integral condition imposed on P1P gives Z P1P (x) dx 0= Ω Z a r2 = 4π c2 + r2 dr 2a 0 3 c2 a a5 = 4π + 3 10a by a change to spherical coordinates. Then c2 = − 3 a 10 and so 3 | x |2 − a. (3.69) 2a 10 We have reached the analogous step to equation 3.51, a two term expansion for the splitting probability: P1P ( x) = " 2 a1 | x| 3 λ0 P (x) = − a + 4πG (x; x1 ) + ǫa1 − N ā + λ0 a2 N ā + λ0 a2 2a 10 (3.70) # N 1 2 a j G x; x j + ǫχ R + O ǫ . − N ā + λ0 a2 ∑ j=1 The derivation of χ R has been left for Appendix A. See Figures 3.17 and 3.18 for some example plots of P . 69 3.3. Multiple targets with asymptotically null volume −3 3 −3 x 10 7 6 2 5 1.5 4 P P 2.5 x 10 1 3 0.5 2 0 1 −0.5 0 0.2 0.4 0.6 0.8 0 0 1 0.2 r 0.4 0.6 0.8 1 r Figure 3.17: Splitting probability from equations 3.70, A.11, and A.19, with ǫ = 0.001 and λ0 = 0.5. P is plotted as a function of r along the radius shown in red (left) or black (right) in Figure 3.8 (the angles which define these radii were chosen randomly). −3 3.5 −4 x 10 10 x 10 9 3 8 7 P P 2.5 2 6 5 1.5 4 1 0.5 0 3 1.5708 3.1416 Azimuthal angle 4.7124 6.2832 2 0 1.5708 3.1416 4.7124 6.2832 Azimuthal angle Figure 3.18: Splitting probability from equations 3.70, A.11, and A.19, with ǫ = 0.001 and λ0 = 0.5. P is plotted as a function of the azimuthal angle along the circle shown in red (left) or black (right) in Figure 3.8 (the radius and the polar angle which define these circles were chosen randomly). 70 3.4. Space dependant diffusivity An analogous derivation yields the MFPT to any absorbing boundary, that is the solution to equation 3.52 changing the Neumann boundary condition by the Robin boundary condition we proposed here ∂n T G = ǫλ0 T G . Because this general MFPT is not what we are interested in, I did not include that calculation. 3.4 Space dependant diffusivity We now consider the problem of space dependant diffusion. As it has been said, lymph nodes have a highly organized architecture, and so it makes sense to think that T cells do not perform free diffusion, but that their diffusion coefficient depends on their location inside the LN. Moreover, it has been proposed [1] that the fibroblastic reticular network might be of aid in the search of a T cell for its cognate APC. In the FRN model dendritic cells are located on the network and T cells preferentially crawl along it with their velocities being higher on the FRN [12]. Figure 3.19: Visual representation of a two dimensional version of our FRN model (the FRN fibers shown in color). The diffusion coefficient on the x (respectively y) direction increases when the cell is on the blue (respectively red) layers. The purple squares mark the intersections between horizontal and vertical fibers of the FRN. On those regions diffusion is higher in both the x and y direction. White squares are areas of slower diffusion off the FRN. This image was adapted from [9]. To mimic the influence of the FRN in the time taken for a T cell to find its corresponding APC, we created a mathematical model where the lymph node is layered in all three dimensions. Diffusion on the x direction (respectively y, z directions) is higher when the cell is on the strips with y and z constant (respectively x and z constant, and x and y constant). (See Figure 3.19 for a visual representation of this model.) 71 3.4. Space dependant diffusivity We write the already familiar equation for the MFPT when there are several absorbing sections on the boundary, in the following way Lǫ T ǫ P = P (x) , (3.71) with P satisfying the corresponding equation from the previous sections. Here ǫ represents the width of the layers. Before we had L = − D ∇2 . This time each entry of the Laplacian will be multiplied by a different diffusion coefficient, which will also depend on x. We want to find T ǫ in the homogenization limit ǫ → 0, i.e. when the strips are so thin that every region, no matter how small, intersects both types of layers. Notice that Lǫ can be written as Lǫ = − A ( x/ǫ) : D 2 where D 2 is the matrix of second partial derivatives, i.e. with i, j entry equal to D2 i,j = ∂2 , ∂xi x j and the colon operator acts in the following way M : N = tr M T N = 3 ∑ Mi,j Ni,j i,j=1 for any matrices M and N. To start with a simpler example, we take a step back to the case where there is only one target. Then the MFPT (to any absorbing part of the boundary) satisfies Lǫ T ǫ = 1 (3.72) [16], where Lǫ = − 3 ∑ i,j=1 Ai,j x ∂2 , ǫ ∂xi x j (3.73) A being the matrix of diffusion coefficients and ǫ the width of the strips. According to [9], in the limit when ǫ → 0 solutions of the general equation Lǫ uǫ = f ( x) converge to solutions of the averaged equation − Ā : D 2 u = f ( x) (3.74) Here the averaged coefficient matrix Ā Ris given by Ā = hρ∞ Ai where the angle brackets denote the average hvi = v, and ρ is the invariant distribution of L, that is the solution to L∗ ρ∞ = 0, hρ∞ i = 1. 72 3.4. Space dependant diffusivity Notice that the matrix A will be periodic since we want the diffusion coefficient to depend only on whether the particle is on or off the color strips and not of on which particular strip it is currently found. For the previous result to hold we need A to have period 1, so we will enforce this when defining the diffusion coefficents on the layered pattern. Thus to compute the average hvi we need only integrate over [0, 1] . To define A, we first consider a symmetric function g ( x) such that g takes higher values on the interval [0, 0.5] , which would correspond to one of the FRN fibers, and lower values on the interval [0.5, 1] . One very simple example is shown in Figure 3.20. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 3.20: Plot of g(z) = (sin (2πz) + 1.1) /2.1. The red line is the border line between layers like in Figure 3.19, the left side of the red line corresponding to a color strip. In terms of such a function g we could define the diffusion coefficient matrix to be g(x2 )g(x3 − 0.5) 0 0 , A ( x) = D 0 g(x1 − 0.5)g(x3 − 0.5) 0 0 0 g(x1 − 0.5)g(x2 ) where x = ( x1 , x2 , x3 ) and D is a basal diffusion coefficient. Clearly A is diagonal and has separable entries, that is it can be expressed as a product of diagonal matrices each depending on only one xi : n A (x) = D ∏ A j x j . j=1 73 3.4. Space dependant diffusivity In the particular case we are considering 1 0 0 1 0 g(x1 − 0.5) 0 A ( x1 ) = 0 0 g(x1 − 0.5) g(x2 ) 0 0 1 0 , A 2 ( x2 ) = 0 0 0 g(x2 ) g(x3 − 0.5) 0 0 0 g(x3 − 0.5) 0 A 3 ( x3 ) = 0 0 1 , . Hence according to [9] the matrix of coefficients for the averaged equation has the simplified form Ā = * 1 j D ∏ nj=1 A jj x j = h Ai . A +* 1 j D ∏nj=1 A jj x j +−1 (3.75) If we use g as in Figure 3.20 then h= Z 1 0 g (z) dz = Z 1 0 g (z − 0.5) dz = 1.1 2.1 so that the averaged equation ends up being − Dh2 ∇2 u = f ( x) (3.76) Notice that the case will be the same whenever the functions that describe the behaviour along each fiber integrate to the same constant, provided that A is diagonal and has separable entries. Another interesting case to consider would be to only have these fast fibers along two of the three directions. Then we would have a different diffusion coefficient along one direction, but all coefficients would still be constant. More complicated dynamics arise if diffusion changes from its basal state in every direction regardless of the orientation of the fiber the cell is at. Although these alternative model is plausible, the simple model presented here already incorporates the influence of the FRN so that there is no need to complicate it. Indeed, by increasing the diffusion coefficient of the same direction as the orientation of the fiber the cell is at, it would appear that it is crawling along the fiber. 74 3.5. Discussion and future work 3.5 Discussion and future work Notice how the MFPT increases 10 fold from a Robin boundary (with λ = 0.5) and 100 fold from a Dirichlet boundary, to a Neumann boundary (see Figure 3.21). Also, in the case where the cell is unable to leave the LN, the MFPT increases very rapidly as the starting position of the T cell moves away from the target. In the other two cases, that is when the cell has some chance to escape from the LN, the increase is less steep. Recall that, in reality, T cells remain in the LN only for some time and if they do not get activated they recirculate through the body and return to the LN later. (a) Neumann (b) Dirichlet (c) Robin Figure 3.21: Mean first passage time with different conditions on the boundary of the LN ∂Ω. Figures 3.1, 3.3 and 3.6 are here shown side by side for ease of comparison. For the case where we considered multiple asymptotically small targets, notice how changing the Neumann for a Robin BC introduces a difference in the splitting probability already at first order. This is so despite the fact that the exit flux out of the LN imposed by the Robin condition is of order ǫ. (Compare equations 3.51 and 3.70.) Indeed, P0 for the Robin condition has an extra λ0 a2 in the denominator so that the probability of hitting Ω1 before the process ends (by absorption by either the other targets or the LN boundary) is smaller in the Robin case. To our knowledge, the splitting probability had never been computed for the multiple asymptotically small tagets with a Robin BC on the outer sphere. That is an important contribution from this thesis work, even if the case we considered is a Robin BC of order ǫ. Moreover, further extensions of this condition can be adapted from the results in [8]. For example, conditions could be derived for a model where the moving particle waits for some time on the boundary before deciding whether to exit or to re75 3.5. Discussion and future work turn to the inside. Such a model might be relevant for other FPT biological problems, if not for the T cell activation problem here proposed. Most results in this chapter can be improved to resemble more closely the dynamics of T cell activation. Many calculations were left in a sense incomplete, for example on Section 3.2.3 I did not compute the variance. Just the calculation of the second moment is very lengthy, complicated and error prone, which is why I did not continue to the variance. In the previous section one could also find explicit solutions for the case when the averaged diffusion coefficients are not equal in all directions. Another interesting variation of the problems presented here would be to include more information about the antigen presenting cells. For example, since dendritic cells are a major class of APCs one could study how having long dendrites changes the time it takes for a T cell to find its corresponding DC. Also all the models presented could be used with real data to estimate contact rates and other kinetic parameters. As I said in the introduction, this mean first passage time approach to the immunological problem was inspired by the Graw-Regoes papers [11, 12]. Hence it is of high interest to compare our results with those presented in the paper. However, this work is not yet in a stage where such comparison can be done. On the 2009 paper, the authors studied the influence of target localization on the killing rates of cytotoxic T cells. Although we considered a different process, namely localization of the activating APC, the problems are very similar and many of the results here presented could be adapted to the other problem. The problem is, however, that our model ends when the target is found. To evaluate the simulated results from the paper under our theoretical framework, we would need to consider a model in which the bond remains for a certain time and after breakage the T cell goes back to the searching stage. I mentioned before that, given the serial engagement hypothesis, this model would also be interesting under the T cell activation problem. Following the FRN model, one could also consider the MFPT to a specific node on a graph (the FRN) when the T cell is only moving along the edges of the graph. This has been studied in the past via simulations [7]. However, a MFPT calculation on a random graph (or other type of graph which might resemble the FRN) is very complicated. In summary, there is a lot more work that can be done and many directions in which this research can be extended. Moreover, I believe the work here presented is significative for it sets the theoretical framework which can not only be used for the immunological problem I proposed, but for 76 3.5. Discussion and future work many other biological problems. 77 Bibliography [1] Marc Bajénoff, Nicolas Glaichenhaus, and Ronald N. Germain. Fibroblastic reticular cells guide T lymphocyte entry and migration within the splenic T cell zone. Journal of Immunology, 181(6), September 2008. [2] Rainer Burkard, Mauro Dell’Amico, and Silvano Martello. Assignment Problems. SIAM, 2009. [3] Claire Chevalier, Olivier Bénichou, Bob Meyer, and Raphaël Voituriez. First-passage quantities of brownian motion in a bounded domain with multiple targets: a unified approach. Journal of Physics A, Mathematical and Theoretical, 44(2), December 2011. [4] Alexei F. Cheviakov and Michael J. Ward. Optimizing the principal eigenvalue of the laplacian in a sphere with interior traps. Mathematical and Computer Modeling, 53(7-8), March 2010. [5] Toby Collins. ordfilt3: Performs 3D order-statistic filtering on 3D volumetric data. From http://www.mathworks.com/matlabcentral/ fileexchange/22044-ordfilt3, 2008. [6] Winfried Denk, James H. Strickler, and Watt W. Webb. Two-photon laser scanning fluorescence microscopy. Science, 248(4951), April 1990. [7] Graham M. Donovan and Grant Lythe. T cell movement on the reticular network. Journal of Theoritical Biology, 295, February 2012. [8] William Feller. Diffusion processes in one dimension. Transactions of the American Mathematical Society, 77(1), May 1954. [9] Brittany D. Froese and Adam M. Oberman. Numerical averaging of non-divergence structure elliptic operators. Communications in Mathematical Sciences, 7(4), August 2009. [10] Izrail’ M. Gradshteyn and Iosif M. Ryzhik. Table of Integrals, Series, and Products. Academic Press, 1980. 78 Bibliography [11] Frederik Graw and Roland R. Regoes. Investigating CTL mediated killing with a 3D cellular automaton. PLoS Computational Biology, 5(8), August 2009. [12] Frederik Graw and Roland R. Regoes. Influence of the fibroblastic reticular network on cell-cell interactions in lymphoid organs. PLoS Computational Biology, 8(3), March 2012. [13] Khuloud Jaqaman, Dinah Loerke, Marcel Mettlen, Hirotaka Kuwata, Sergio Grinstein, Sandra L. Schmid, and Gaudenz Danuser. Robust single-particle tracking in live-cell time-lapse sequences. Nature Methods, 5(8), 2008. Available at http://lccb.hms.harvard.edu/ software.html. [14] Yoshitsugu Kawashima, Makoto Sugimura, Yann-Ching Hwang, and Norio Kudo. The lymph system in mice. The Japanese Journal of Veterinary Research, 12(4), July 1964. [15] Theodore Kolokolnikov, Michèle S. Titcombe, and Michael J. Ward. Optimizing the fundamental Neumann eigenvalue for the Laplacian in a domain with small traps. European Journal of Applied Mathematics, 16(2), April 2005. [16] Hannah W. McKenzie, Mark A. Lewis, and Evelyn H. Merrill. First passage time analysis of animal movement insights into the functional response. Bulletin of Mathematical Biology, 71(1), September 2009. [17] MicroscopyU. Fundamentals and applications in multiphoton excitation microscopy. From http://www.microscopyu.com/articles/ fluorescence/multiphoton/multiphotonintro.html. [18] Jennifer S. Morrison. Deciphering multi-state mobility within single particle trajectories of proteins on the plasma membrane. Master’s thesis, The University of British Columbia, August 2010. [19] Kenneth Murphy. Janeway’s Immunobiology. Garland Science, 8th edition, 2012. [20] National Institutes of Health. Image J. From http://rsbweb.nih. gov/ij/. [21] Jean-Christophe Olivo-Marin. Extraction of spots in biological images using multiscale products. Pattern Recognition, 35(9), 2002. 79 [22] Sidney Redner. A Guide to First passage processes. Cambridge University Press, 2001. [23] Hannes Risken. The Fokker-Planck Equation. Methods of Solution and Applications. Springer-Verlag, second edition, 1989. [24] Ulrich H. von Andrian and Thorsten R. Mempel. Homing and cellular traffic in lymph nodes. Nature Reviews Immunology, 3:867–878, November 2003. [25] Michael J. Ward. Math 551: Homework 2. From http://www.math. ubc.ca/~ward/teaching/m551/hw2_11sol.pdf. [26] Michael J. Ward. Notes on strong localized perturbation theory. From http://www.math.ubc.ca/~ward/teaching/m551/m550_hole.pdf. [27] Carla Wofsy, Daniel Coombs, and Byron Goldstein. Calculations show substantial serial engagement of T cell receptors. Biophysical Journal, 80(2), February 2001. 80 Appendix A Additional derivations A.1 From the Neumann splitting probability, equation 3.51 A.1.1 The Neumann Green’s function for the sphere We now calculate the Neumann Green’s function G ( x, ξ ) satisfying equation 3.49 when Ω is a sphere of radius a. First let us look at the singular part of G, that is we want to solve ∇2 u (x) = −δ (x − ξ) , in R3 . (A.1) Take a small sphere Ωǫ of radius ǫ about ξ. Then in a neighborhood of ξ we let r = | x − ξ | and then 0 = ∇2 u 2 = urr + ur , for r > 0, r β so that u (r) = r . Now, since Z Ωǫ ∇2 u ( x) dx = − Z Ωǫ δ ( x − ξ ) dx = −1, then by the divergence theorem −1 = Z Ωǫ ∇u ·~ndS = 4π r2 ur (r) |r=ǫ , (A.2) again by a change to spherical coordinates as in equation 3.45. Since ur = − βr−2 then equation A.2 gives ǫ2 β 4π − 2 ǫ = −1 81 A.1. From the Neumann splitting probability, equation 3.51 which yields β= Thus u (x) = 1 . 4π 1 , 4π | x − ξ | and so we can write G (x; ξ) ≈ 1 + R (x; ξ) , as x → ξ. 4π | x − ξ | (A.3) To find R I followed the Appendix A from reference [15]. We first decompose G as follows G ( x; ξ ) = u ( x) + | x |2 + | ξ |2 1 + φ ( x; ξ ) , 6| Ω | 4π (A.4) with φ an unknown function we want to compute. Notice that in that case 2 2∂ ∂ 1 1 2 2 2 2 2 ∇ φ (x; ξ) + ∇ G (x; ξ) = ∇ u (x) + | x | + | ξ | + 2 6|Ω| ∂r r∂ 4π 1 2 1 + ∇ φ (x; ξ) . = −δ (x − ξ) + |Ω| 4π Because Ω is a sphere, |Ω| = 4πa3 /3. Then 1 1 | x| ∂n G (x; ξ) = ∂n + 3 + ∂n φ . 4π |x − ξ | a Substituting this onto the differential equation which defines G, equation 3.49, we find that φ must satisfy ∇2 φ = 0, x ∈ Ω, 1 1 ∂n φ = − 2 − ∂n , x ∈ ∂Ω. a |x − ξ | (A.5) To solve A.5 we choose a coordinate system so that the source point x = ξ is on the positive z axis. Then, since ∇2 φ = 0 and φ is axisymmetric, φ admits the series expansion ∞ | x||ξ | n , (A.6) φ ( x; ξ ) = ∑ Bn Pn (cos θ ) a2 n=0 82 A.1. From the Neumann splitting probability, equation 3.51 where Pn is the Legendre polynomial of integer n and the Bn for n = 0, 1, . . . are coefficients to be determined. To simplify the problem, we have enforced that φ ( x; ξ ) = φ (ξ; x) so that G ( x; ξ ) = G (ξ; x), an assumption that agrees with the interpretation of P . From A.6 we see that ∞ ∂n φ|∂Ω = nBn ∑ an+1 Pn (cos θ ) |ξ |n . (A.7) n=0 Now, from the generating function for Legendre polynomials, √ 1 1 − 2tz + t2 ∞ = ∑ tk Pk (z) , (A.8) k=0 we obtain the following equation, by setting z = cos θ and t = |ξ |/| x|: 1 1 = p 2 |x − ξ | |ξ | − 2| x||ξ | cos θ + | x|2 1 √ = 2 | x| t − 2tz + 1 ∞ |ξ |n =∑ P (cos θ) . n+1 n n=0 | x | By differentiating this expression with respect to r = | x| we get ∞ 1 (n + 1) |ξ |n ∂r Pn (cos θ) . |r=a = − ∑ |x − ξ | an+2 n=0 (A.9) Upon substituting A.7 and A.9 into the boundary condition given in A.5, we obtain ∞ 1 n+1 |ξ |n ∑ nBn − a Pn (cos θ ) an+1 = − a2 . n=0 Since B0 is arbitrary, we can choose it to be equal to 1/a for convenience. Then the other coefficients must satisfy Bn = 1 1 + , for n ≥ 1. a na (A.10) Back to the series expansion of φ, equation A.6, and replacing on it the previous expressions for the Bn ’s, we find that | x||ξ | n 1 ∞ 1 | x||ξ | n 1 ∞ + . θ) P φ (x; ξ) = ∑ Pn (cos θ) ∑ n n (cos a n=0 a2 a n=1 a2 83 A.1. From the Neumann splitting probability, equation 3.51 By using the generating function A.8 again we find that the first infinite sum reduces to 1 ∞ Pn (cos θ) a n∑ =0 | x||ξ | a2 n 4 − 12 a a a2 2 = − 2| ξ | cos θ + |ξ | | x | | x |2 | x| a , = | x |r ′ where x′ = xa2 /| x|2 is the image point to x outside the sphere and r′ = | x ′ − ξ |. To calculate the second infinite sum let us define the following function: ∞ I ( β) = 1 ∑ n Pn (cos θ ) βn . n=1 In terms of it, the second infinite sum is equal to a−1 I | x||ξ |/a2 . Upon diferentiating and then using the generating function A.8 we get ∞ I′ β = = ∑ Pn (cos θ) βn−1 n=1 1 β 1 = β 1 = β ∞ ∑ Pn (cos θ) βn n=1 ∞ ∑ βn Pn (cos θ) − P0 (cos θ) n=0 ! ! 1 p −1 , β2 − 2β cos θ + 1 since P0 (cos θ ) = 1. We can integrate the previous equation [see 10, page 95] to obtain an expression for I. Since I (0) = 0, Z β1 1 √ − 1 ds I β = 0 s s2 − 2s cos θ + 1 ! 2 p = log . 1 − β cos θ + 1 + β2 − 2β cos θ We are finally in a position to obtain a final expression for φ and so also for G. a 2a2 1 φ (x; ξ) = + log , | x |r ′ a a2 − | x||ξ | cos θ + | x|r′ 84 A.1. From the Neumann splitting probability, equation 3.51 | x |2 + | ξ |2 a 1 + + 4π | x − ξ | 6| Ω | 4π | x|r′ (A.11) 2a2 1 + log + B, 4πa a2 − | x||ξ | cos θ + | x|r′ R where B is a constant which willRbe chosen to satisfy Ω G = 0. Now we want to show that Ω G ( x; ξ ) is independent R of ξ and so that we can compute B by looking at the simpler equation Ω G ( x; 0) = 0. We have: Z 1 ′ 2 + δ (x − ξ) dx G x; ξ 0= ∇ G (x; ξ) − |Ω| Ω Z Z 1 G x; ξ ′ ∇2 G (x; ξ) dx − = G x; ξ ′ + G ξ; ξ ′ |Ω| Ω Z ZΩ ′ ∇ G (x; ξ) · ∇ G x; ξ ′ dx ∇ · G x; ξ ∇ G (x; ξ) dx − = G (x; ξ) = Ω Ω 1 − |Ω| Z Ω G x; ξ ′ + G ξ; ξ ′ . Then the divergence theorem gives Z Z ′ ∇ · G x; ξ ∇ G (x; ξ) dx = Ω ∂Ω = 0, G x; ξ ′ ∇ G (x; ξ) · ~ndx by the boundary condition on G, i.e. ∂n G = 0 on ∂Ω. Hence Z Z 1 ∇ G (x; ξ) · ∇ G x; ξ ′ dx. G x; ξ ′ = G ξ; ξ ′ − |Ω| Ω Ω R Note that the RHS is symmetric in ξ and ξ ′ . It follows that Ω G ( x; ξ ) dx = R ′ Ω G ( x; ξ ) dx. With a change to spherical coordinates, like in equation 3.45, we compute the desired integral: 0= Z Ω 1 = 4π G (x; 0) dx Z 2π Z π Z ∞ 0 1 + 6| Ω | Z a 0 0 1 r sin θdrdθdϕ + 4πa Z 2π Z π Z ∞ 0 0 a 0 0 0 0 r2 sin θdrdθdϕ 4 r sin θdrdθdϕ + B|Ω| Z a 1 4π = rdr + r2 dr + a 6 |Ω| 0 0 1 1 1 4πa3 + + B. + = a2 2 3 10 3 Z Z 2π Z π Z ∞ 0 r4 dr + B|Ω| 85 A.1. From the Neumann splitting probability, equation 3.51 Here we have used the fact that when ξ = 0, r′ = | x′ | = a2 /| x|. Thus B=− 7 . 10πa (A.12) A.1.2 Derivation of χ1 In this section we compute the constant from the order ǫ function, P1 , from the splitting probability with multiple targets and Neumann outer boundary. Let us first write P1 ( x) as x → x j using what we found on the previous section (equation A.3). P1 (x) ≈ N (1 − P0 ) a1 + 4π (1 − P0 ) a1 R1,1 − 4π P0 ∑ ai G1,i + χ1 , as x → x1 , | x − x1 | i=2 P1 (x) ≈ − N P0 a j − 4π P0 a j R j,j + 4πa1 Gj,1 − 4π P0 ∑ ai Gj,i + χ1 , |x − xj | i=1,i6= j as x → x j , j = 2, . . . , N, where Gi,j ≡ G xi ; x j and Ri,j ≡ R xi ; x j . We reduce the previous equations by introducing the notation ! N B1 = 4πa1 R1,1 − 4π P0 a1 R1,1 + ∑ ai G1,i , i=2 N Bj = 4πa1 Gj,1 − 4π P0 Thus a j R j,j + ∑ ai Gj,i i=1,i6= j ( 1 − P 0 ) a1 + B1 + χ1 , | x − x1 | P1 ( x ) ≈ P0 a j + B j + χ1 , − | x − x1 | ! , j = 2, . . . , N. as x → x1 . (A.13) as x → x j From the matching conditions 3.41 and 3.42 we get that w1 must satisfy w1 ≈ Bj + χ1 as |y| → ∞. On the other hand, equation 3.40 has solution aj , w1 ( y) = c 1 − |y| where c is a constant which is readily found to be c = Bj + χ1 , since w1 → c as |y| → ∞. 86 A.1. From the Neumann splitting probability, equation 3.51 Going back again to the matching condition 3.41, it yields a j B j + χ1 P2 ≈ − , as x → x j , j = 1, . . . , N. |x − xj | By an argument like that led by equations 3.44 and 3.45 we find that P2 should satisfy N ∇2 P2 = 4π ∑ a j Bj + χ1 δ x − x j , x ∈ Ω, (A.14) j=1 ∂n P2 = 0, x ∈ ∂Ω. The divergence theorem yields a solvability condition for P2 , just as it did for P1 on equation 3.47, namely N ∑ aj j=1 Bj + χ1 = 0. Hence χ1 = − 1 N ā N ∑ a j Bj . j=1 The expressions for Bj can be substituted back and so χ1 can be written in matrix form as 4πa1 T 1 T χ1 = − G a − a Ga , (A.15) N ā N ā 1 where a = ( a1 , . . . , a N ) T and R1,1 G1,2 ··· G1,N .. . . .. .. G2,1 . G≡ . . . .. .. .. GN −1,N GN,1 · · · GN,N −1 R N,N . Notice that with equations A.11 and A.12 we can compute all the entries of this matrix. In particular a2 1 a 7 | x i |2 + + log . (A.16) − Ri,i = 4πa3 4π a2 − | xi |2 4πa a2 − | x i |2 10πa 87 A.2. From the Robin splitting probability, equation 3.70 A.2 From the Robin splitting probability, equation 3.70 A.2.1 Derivation of χ R In this section we compute the constant from the order ǫ function, P1 , from the splitting probability with multiple targets and Robin outer boundary. First we look at the limit as x → x j of P1 ( x). Like in Section A.1.2 we write ( 1 − P 0 ) a1 + B1 + χ R , as x → x1 | x − x1 | P1 ( x ) ≈ (A.17) P0 a j + Bj + χ R , as x → x j − | x − x1 | with N a1 R1,1 + ∑ ai G1,i B1 = −λ0 P0 P1P (x1 ) + 4πa1 R1,1 − 4π P0 i=2 ! , ! , j = 2, . . . , N. N Bj = −λ0 P0 P1P x j + 4πa1 Gj,1 − 4π P0 a j R j,j + ∑ i=1,i6= j ai Gj,i Following an argument entirely analogous to that on A.1.2 we find that P2 must satisfy N ∇2 P2 = 4π ∑ a j Bj + χ R δ x − x j , x ∈ Ω, (A.18) j=1 ∂n P2 = −λ0 P1 , x ∈ ∂Ω. Upon using the divergence theorem on the previous equation we obtain a solvability condition for P2 , namely N 4π Z B + χ a ∑ j j R = ∇2 P2 dV j =1 Ω = − λ0 Z = λ20 P0 Z∂Ω P1 dS ∂Ω P1P dS N − λ0 4π ∑ δj1 − P0 a j P1P x j − λ0 χ R |∂Ω| j=1 88 A.2. From the Robin splitting probability, equation 3.70 by equation 3.68. Now notice that P1P is constant on ∂Ω since it is a function of only | x|. Its value is a/2 − 3a/10 = a/5. Hence Z ∂Ω P1P dS = a|∂Ω| 4πa3 = . 5 5 Substituting back the expressions for Bj we find that ∑ N j=1 a j B j can be written with matrix notation as N h i N T T a B = − λ P x + 4π a G a a P − P a G a . o 0 ∑ j 1P 0 1 j ∑ j j j =1 1 j=1 Recall that P1P ( x) = 3 | x |2 − a. 2a 10 Therefore 1 χR = 4π N ā + λ0 a2 λ20 P0 Z ∂Ω N P1P dS − λ0 4π ∑ δj1 − P0 a j P1P x j j=1 N − 4π ∑ a j Bj j=1 = 1 N ā + λ0 a2 N 1 2 λ0 P0 a3 − λ0 a1 P1P (x1 ) + λ0 P0 ∑ a j P1P x j 5 j=1 N h T + λo P0 ∑ a j P1P x j − 4π a1 G a j=1 = 1 N ā + λ0 a2 ! λ20 P0 a3 λ0 a1 | x1 |2 3λ0 a1 a λ0 P0 − + + 5 2a 10 a N 1 T − P0 a G a i ∑ a j | x j |2 j=1 i 3λ0 P0 aN ā − − 4π a1 G T a − P0 aT G a 5 1 h ! ! . (A.19) 89 Appendix B Programs In this appendix I present the most relevant codes written by me which I used for Chapter 2. 90 B.1 spotDetector3D Extension to three dimensions of the program spotDetector, algorithm from [21], coded by François Aguet, which can be downloaded from http://lccb.hms.harvard.edu/software.html. B.1. spotDetector3D 91 function movieInfo=spotDetector3D(path,Z,T,choice,dthreshold) % Adapted for 3D from SpotDetector, http://lccb.hms.harvard.edu/doc/spotDetector_101410.zip % % Reads a denoised 3D video in the form of a collection of 2D images % Detects centers of 3D cells in each snapshot % Centers are determined in the following way: % - First, find all 3D connected regions (considering only nearest % neighbors). % - Get the centers of mass of these regions. % - Find local maxima (using a cubic window of side 7). % - For each region, retain % - the first local maximum if choice=0 % - the center of mass if choice=1 % and the secondary local maxima that are farther away from both the % center of mass and the highest local maximum than dthreshold. % % INPUT path is the path to the folder where the images are contained % the code assumes that such folder follows the pattern created by % imstotiff.m, i.e. % - such folder has the name of the data set, with which also the images % names start % - inside this folder, there’s another one named "processed" where all the denoised images are contained - the images are named name_Tddd_Zdd.tif Images are also assumed to be squared Z is the number of planes from a 3D "snapshot" T is the total number of time points imaged OUTPUT is movieInfo, appropiately formatted to be input to u-track For a movie with N frames, movieInfo is a structure array with N entries. Every entry has the fields xCoord, yCoord, zCoord (if 3D) and amp. If there are M objects in frame i, each one of these fields in moveiInfo(i) will be an Mx2 array, where the first column is the value (e.g. x-coordinate in xCoord and amplitude in amp) and the second column is all zero. In SpotDetector3D, flag is unused. Uncomment the block starting with "if flag==0", to split also regions where cells "look like" 3D 8’s. 92 if nargin<5 dthreshold = 5; end if nargin<4 choice=0; elseif isempty(choice) choice=0; end B.1. spotDetector3D % % % % % % % % % % % % % % % % % movieInfo=repmat(struct(’xCoord’,[],’yCoord’,[],’zCoord’,[],’amp’,[]),T,1); [~,base]=fileparts(path); L=size(imread([path,filesep,’processed’,filesep,base,’_T001_Z01.tif’]),1); img3d=zeros(L,L,Z); %img3d will contain the 3D image ny=L; nx=L; nz=Z; %% 3D process connected components mask3d=img3d>0; localMax = locmax3d(img3d, 7); [labels, nComp] = bwlabeln(mask3d, 6); 93 area = zeros(nComp, 1); B.1. spotDetector3D for t=1:T % save 2D images into 3D compund tname=[’_T’,num2str(t,’%3.3d’)]; for z=1:Z sufix=[’_Z’,num2str(z,’%2.2d’),’.tif’]; name=[path,filesep,’processed’,filesep,base,tname,sufix]; img3d(:,:,z)=imread(name); end totalInt = zeros(nComp, 1); nMaxima = zeros(nComp, 1); xmax = zeros(nComp, 1); ymax = zeros(nComp, 1); zmax = zeros(nComp, 1); xcom = zeros(nComp, 1); ycom = zeros(nComp, 1); zcom = zeros(nComp, 1); labelVect = zeros(nComp, 1); % Compute area and center of mass for each component stats = regionprops(labels, img3d, ’Area’, ’WeightedCentroid’, ’PixelIdxList’); % Component labels of local maxima maxLabels = labels .* (labels & localMax>0); maxCoords(1:nComp) = struct(’PixelIdxList’, []); mc = regionprops(maxLabels, ’PixelIdxList’); maxCoords(1:length(mc)) = deal(mc); B.1. spotDetector3D xmax2 = cell(nComp, 1); ymax2 = cell(nComp, 1); zmax2 = cell(nComp, 1); area2 = cell(nComp, 1); totalInt2 = cell(nComp, 1); labelVect2 = cell(nComp, 1); 94 for n = 1:nComp [yi,xi, zi] = ind2sub([ny nx nz], stats(n).PixelIdxList); [ym,xm,zm] = ind2sub([ny nx nz], maxCoords(n).PixelIdxList); area(n) = stats(n).Area; com = stats(n).WeightedCentroid; xcom(n) = com(1); ycom(n) = com(2); zcom(n) = com(3); 95 nMaxima(n) = length(xm); if nMaxima(n)==1 xmax(n) = xm; ymax(n) = ym; zmax(n) = zm; nMaxima(n) = 1; labelVect(n) = labels(ym,xm,zm); elseif nMaxima(n)==0 % no maximum was detected for this cluster maxValueIdx = find(values == max(values)); xmax(n) = xi(maxValueIdx(1)); ymax(n) = yi(maxValueIdx(1)); zmax(n) = zi(maxValueIdx(1)); nMaxima(n) = 1; B.1. spotDetector3D values = img3d(stats(n).PixelIdxList); totalInt(n) = sum(values); labelVect(n) = labels(ymax(n), xmax(n), zmax(n)); else % resolve multiple maxima cases maxValues = localMax(sub2ind(size(localMax), ym, xm, zm)); % highest local max maxIdx = find(maxValues == max(maxValues)); xmax(n) = xm(maxIdx(1)); ymax(n) = ym(maxIdx(1)); zmax(n) = zm(maxIdx(1)); labelVect(n) = labels(ymax(n), xmax(n), zmax(n)); % compute distance of secondary maxima to primary dist2max = sqrt((xmax(n)-xm).^2 + (ymax(n)-ym).^2 + (zmax(n)-zm).^2); dist2com = sqrt((xcom(n)-xm).^2 + (ycom(n)-ym).^2 + (zcom(n)-zm).^2); mindist = min(dist2max,dist2com); 96 % retain secondary maxima where mindist > threshold idx2 = find(mindist > dthreshold); if ~isempty(idx2) xmax2{n} = xm(idx2); ymax2{n} = ym(idx2); zmax2{n} = zm(idx2); nSecMax = length(idx2); B.1. spotDetector3D % remove highest max from list xm(maxIdx(1)) = []; ym(maxIdx(1)) = []; zm(maxIdx(1)) = []; nMaxima(n) = nSecMax+1; % split area area2{n} = area(n)*ones(nSecMax,1)/nMaxima(n); area(n) = area(n)/nMaxima(n); labelVect2{n} = labels(sub2ind(size(labels), ymax2{n}, xmax2{n}, zmax2{n})); end end end xmax2 = vertcat(xmax2{:}); ymax2 = vertcat(ymax2{:}); zmax2 = vertcat(zmax2{:}); totalInt2 = vertcat(totalInt2{:}); area2 = vertcat(area2{:}); labelVect2 = vertcat(labelVect2{:}); 97 % assign results frameInfo.xmax = frameInfo.ymax = frameInfo.zmax = frameInfo.xcom = to output structure [xmax; xmax2(:)]; [ymax; ymax2(:)]; [zmax; zmax2(:)]; [xcom; xmax2(:)]; B.1. spotDetector3D %intensity values totalInt2{n} = totalInt(n)*ones(nSecMax,1)/nMaxima(n); totalInt(n) = totalInt(n)/nMaxima(n); frameInfo.ycom = [ycom; ymax2(:)]; frameInfo.zcom = [zcom; zmax2(:)]; frameInfo.totalInt = [totalInt; totalInt2(:)]; frameInfo.area = [area; area2(:)]; frameInfo.nMaxima = nMaxima; % maxima per component frameInfo.labels = [labelVect; labelVect2(:)]; frameInfo.nComp = nComp; 98 movieInfo(t).amp(:,1) = frameInfo.totalInt; if choice==0 movieInfo(t).xCoord(:,1) = frameInfo.xmax; movieInfo(t).yCoord(:,1) = frameInfo.ymax; movieInfo(t).zCoord(:,1) = frameInfo.zmax; elseif choice==1 movieInfo(t).xCoord(:,1) = frameInfo.xcom; movieInfo(t).yCoord(:,1) = frameInfo.ycom; movieInfo(t).zCoord(:,1) = frameInfo.zcom; B.1. spotDetector3D % prepare fields for tracker nObj = length(frameInfo.xmax); movieInfo(t).amp = zeros(nObj,2); movieInfo(t).xCoord = zeros(nObj,2); movieInfo(t).yCoord = zeros(nObj,2); movieInfo(t).zCoord = zeros(nObj,2); end frameInfo.path = []; frameInfo.maskPath = []; end B.2 mySpotDetector Alternative 3D detection algorithm of my own autorship. 99 function movieInfo=mySpotDetector(path,Z,T,bt) % Algorithm by Monica Delgado % Uses ideas from SpotDetector, http://lccb.hms.harvard.edu/doc/spotDetector_101410.zip % % Reads a denoised 3D video in the form of a collection of 2D images % Detects centers of 3D cell-like connected regions in each snapshot % The cell-like characteristic is obtained by splitting regions B.2. mySpotDetector parent=fileparts(pwd); if choice==0 sfx=’max’; elseif choice==1 sfx=’com’; end save([parent,filesep,’Results’,filesep,base,’sd3D_’,sfx,’.mat’],’movieInfo’); that look like two or more sphere-like objects pasted INPUT path is the path to the folder where the images are contained the code assumes that such folder follows the pattern created by imstotiff.m, i.e. - such folder has the name of the data set, with which also the images names start - inside this folder, there’s another one named "processed" where all the denoised images are contained - the images are named name_Tddd_Zdd.tif Images are also assumed to be squared Z is the number of planes from a 3D "snapshot" T is the total number of time points imaged bt is an indicator - if bt=1 then sweeping will be done from top to bottom - if bt=0 then sweeping will be done from bottom to top (default) Results will be saved in folder Results under the main folder Output is movieInfo, appropiately formatted to be input to u-track For a movie with N frames, movieInfo is a structure array with N entries. Every entry has the fields xCoord, yCoord, zCoord (if 3D) and amp. If there are M objects in frame i, each one of these fields in moveiInfo(i) will be an Mx2 array, where the first column is the value (e.g. x-coordinate in xCoord and amplitude in amp) and the second column is all zero. B.2. mySpotDetector % % % % % % % % % % % % % % % % % % % % % % % % % 100 if nargin<4 bt=0; end movieInfo=repmat(struct(’xCoord’,[],’yCoord’,[],’zCoord’,[],’amp’,[]),T,1); [~,base]=fileparts(path); L=size(imread([path,filesep,’processed’,filesep,base,’_T001_Z01.tif’]),1); imagen=zeros(L,L,Z); %imagen will contain the 3D image fprintf(’Starting analysis of time point %d...\n’,t); 101 % begins analysis bimage=imagen>0; %make binary image corresponding to 3d image CC=bwconncomp(bimage,6); % find 3D connected regions R=CC.PixelIdxList; n=CC.NumObjects; B.2. mySpotDetector for t=1:T % save 2D images into 3D compund tname=[’_T’,num2str(t,’%3.3d’)]; for z=1:Z sufix=[’_Z’,num2str(z,’%2.2d’),’.tif’]; name=[path,filesep,’processed’,filesep,base,tname,sufix]; imagen(:,:,z)=imread(name); end stats=regionprops(CC,imagen,’BoundingBox’,’WeightedCentroid’,’Image’); for i=1:n % sweep all connected regions bb=ceil(stats(i).BoundingBox); %the first 3 entries of bb denote the coordinates of the upper %right pixel of the bounding box, the following 3 denote the %width of the bounding box along each dimension no 0]; 0]; 0]; else flag=0; %this flag will turn on if in some plane of the current %3D region there’s more than one 2D connected region 102 s=size(stats(i).Image); img2d=zeros(s); img2d(stats(i).Image==1)=imagen(R{i}); % img2d is the non-binary version of output Image from % regionprops B.2. mySpotDetector if ((bb(4)==1) || (bb(5)==1) ||(bb(6)==1)) % if the region is really of dim less than 3, then there’s % problem determining its center movieInfo(t).xCoord(end+1,:)=[stats(i).WeightedCentroid(1) movieInfo(t).yCoord(end+1,:)=[stats(i).WeightedCentroid(2) movieInfo(t).zCoord(end+1,:)=[stats(i).WeightedCentroid(3) movieInfo(t).amp(end+1,:)=[sum(imagen(R{i})) 0]; 103 % prevPIL=previous Pixel Indexes List, % sweep all subsequent planes included in the 3D region % this algorithm will find the intersections of 2D regions % between planes if bt==1 ss=bb(6)-1; B.2. mySpotDetector % initialize: find connected 2D regions inside the first % plane of the current 3D region if bt==1 cc2d=bwconncomp(stats(i).Image(:,:,end),8); % start on top plane elseif bt==0 cc2d=bwconncomp(stats(i).Image(:,:,1),8); % start on bottom plane end if cc2d.NumObjects>1 flag=1; end prevPIL=cc2d.PixelIdxList; for m=1:length(prevPIL) %for each 2d connected region if bt==1 prevplane{m}=bb(6); % start on top plane elseif bt==0 prevplane{m}=1; % start on bottom plane end end 104 m=1; % new regions counter for j=1:length(prevPIL) flagj=0; % if a given previous 2D region does not have % intersection with regions in the current plane, % don’t intersect it (which would erase it), keep it, as % it is a region which has probably ended in the % previous plane, that’s what this flagj is for % Same way, we must keep the 2D regions in current B.2. mySpotDetector ff=1; inc=-1; elseif bt==0 ss=2; ff=bb(6); inc=1; end for z=ss:inc:ff %find 2D connected regions in each plane of the current %3D region cc2d=bwconncomp(stats(i).Image(:,:,z),8); if cc2d.NumObjects>1 flag=1; end newlist=cc2d.PixelIdxList; % current plane list of pixels in each region flaglist=1:cc2d.NumObjects; % plane which had empty intersections with all previous % regions, flaglist helps here if flagj==0 PIL{m}=prevPIL{j}; plane{m}=prevplane{j}; m=m+1; end end 105 % keep 2D regions in current plane that matched no % previous regions B.2. mySpotDetector %compare each region found so far with %current regions, keep non-empty intersections for k=1:cc2d.NumObjects set=intersect(prevPIL{j},newlist{k}); if ~isempty(set) flagj=1; flaglist=setdiff(flaglist,k); PIL{m}=set; plane{m}=[prevplane{j},z]; m=m+1; end end for k=flaglist PIL{m}=newlist{k}; plane{m}=z; m=m+1; end prevPIL=PIL; prevplane=plane; % % % % % % 106 if flag==0 movieInfo(t).xCoord(end+1,:)=[stats(i).WeightedCentroid(1) 0]; movieInfo(t).yCoord(end+1,:)=[stats(i).WeightedCentroid(2) 0]; movieInfo(t).zCoord(end+1,:)=[stats(i).WeightedCentroid(3) 0]; movieInfo(t).amp(end+1,:)=[sum(imagen(R{i})) 0]; else % convert PIL and plane to 3D linear indeces of Image label=zeros(s); mm=1; for m=1:length(PIL) if ((length(PIL{m})>1) || (length(plane{m})>1)) lidx{mm}=ind2sub2ind(s(1:2),PIL{m},plane{m}); label(lidx{mm})=mm; mm=mm+1; end B.2. mySpotDetector end % end of the current plane of a given 3D region PIL={}; plane={}; lidx={}; end %ends case more than one z-slice end % end of current 3D region end %end of current time image 107 parent=fileparts(pwd); if bt==0 sfx=’bup’; B.2. mySpotDetector % end % find the weighted centroids of the regions found aux=regionprops(label,img2d,’WeightedCentroid’); cntr=zeros(length(aux),3); for m=1:length(aux) cntr(m,:)=aux(m).WeightedCentroid; movieInfo(t).amp(end+1,:)=[sum(img2d(lidx{m})) 0]; end %recover corresponding coordinates in original image movieInfo(t).xCoord(end+1:end+m,:)=[cntr(:,1)+bb(1) zeros(m,1)]; movieInfo(t).yCoord(end+1:end+m,:)=[cntr(:,2)+bb(2) zeros(m,1)]; movieInfo(t).zCoord(end+1:end+m,:)=[cntr(:,3)+bb(3) zeros(m,1)]; end elseif bt==1 sfx=’tdown’; end save([parent,filesep,’Results’,filesep,base,’mysd_’,sfx,’.mat’],’movieInfo’); B.3 detecttr This function uses u-track on the denoised images both for detection of 3D cell centers and for the actual tracking. function H=detecttr(path,Z,T,d) prnt=pwd; [~,name]=fileparts(path); if d==0 mkdir(path,’lessnoise’) end movieInfo=repmat(struct(’xCoord’,[],’yCoord’,[],’zCoord’,[],’amp’,[],’height’,[]),T,1); H=[]; 108 for t=1:T cd([prnt,filesep,’SpotDetector’]) tstr=num2str(t,’%3.3d’); B.3. detecttr if nargin<4 d=0; end B.3. detecttr 109 for z=1:Z if d==0 frame = double(imread([path,filesep,’original’,filesep,name,’_T’,tstr,’_Z’,... num2str(z,’%2.2d’),’.tif’])); denoised=denoise(frame,4,2); imwrite(uint8(denoised-1),[path,filesep,’lessnoise’,filesep,name,’_T’,tstr,... ’_Z’,num2str(z,’%2.2d’),’.tif’],’tiff’); elseif d==1 denoised=imread([path,filesep,’lessnoise’,filesep,name,’_T’,tstr,’_Z’,... num2str(z,’%2.2d’),’.tif’]); end mi(z) = spotDetector(denoised); %H=[H; mi(z).area]; H=[H; mi(z).totalInt./mi(z).area]; end cd([prnt,filesep,’u-track111221’,filesep,’u-track’]) tfdet=scriptTrackGeneral(mi,0,2,1,5,0,1,2); % last parameter is 1 for all LN except z1, in which case it’s 2 for i=1:length(tfdet) if any(isnan(tfdet(i).tracksCoordAmpCG(1:8:end))) print(’Gap found in track %d, t = %d\n’,i,t) print(’*****************************\n’) else movieInfo(t).xCoord=[movieInfo(t).xCoord; mean(tfdet(i).tracksCoordAmpCG(1:8:end)) 0]; movieInfo(t).yCoord=[movieInfo(t).yCoord; mean(tfdet(i).tracksCoordAmpCG(2:8:end)) 0]; zc=sum(tfdet(i).tracksCoordAmpCG(4:8:end).*[tfdet(i).seqOfEvents(1,1):... tfdet(i).seqOfEvents(2,1)])/sum(tfdet(i).tracksCoordAmpCG(4:8:end)); movieInfo(t).zCoord=[movieInfo(t).zCoord; zc 0]; movieInfo(t).height=[movieInfo(t).height; tfdet(i).seqOfEvents(2,1)... -tfdet(i).seqOfEvents(1,1)+1]; movieInfo(t).amp=[movieInfo(t).amp; max([std(tfdet(i).tracksCoordAmpCG(1:8:end)),... std(tfdet(i).tracksCoordAmpCG(2:8:end)),std(tfdet(i).tracksCoordAmpCG(3:8:end))]) 0]; end end end tf=scriptTrackGeneral(movieInfo,0,3,1,2,2,0,10); cd(prnt) %mkdir(’Results’,[name,’trd’]) save([’Results’,filesep,name,’trd’,filesep,name,’trdlnzamp.mat’],’movieInfo’); save([’Results’,filesep,name,’trd’,filesep,’tr’,name,’-trdlnzamp.mat’],’tf’); 110 Function scriptTrackGeneral is called in the following way: B.3. detecttr for i=1:T %movieInfo(i).xCoord(:,1)=movieInfo(i).xCoord(:,1)*530.33/512; %movieInfo(i).xCoord(:,1)=movieInfo(i).xCoord(:,1)*424.27/512; %for LN10 movieInfo(i).xCoord(:,1)=movieInfo(i).xCoord(:,1)*283.4/560; %for z1 movieInfo(i).yCoord(:,1)=movieInfo(i).yCoord(:,1)*283.4/560; movieInfo(i).zCoord(:,1)=2.5*(movieInfo(i).zCoord(:,1)-1); %the scaling factor is 2 for most LN except for z1 in which case it is 2.5 end B.3. detecttr 111 function tracksFinal=scriptTrackGeneral(movieInfo,tosave,probDim,TW,minTL,LM,minSR,maxSR) % Copyright (C) 2011 LCCB % % This file is part of u-track. % % u-track is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % (at your option) any later version. % % u-track is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with u-track. If not, see <http://www.gnu.org/licenses/>. % % %% Detection results: movieInfo % %For a movie with N frames, movieInfo is a structure array with N entries. %Every entry has the fields xCoord, yCoord, zCoord (if 3D) and amp. %If there are M features in frame i, each one of these fields in %moveiInfo(i) will be an Mx2 array, where the first column is the value %(e.g. x-coordinate in xCoord and amplitude in amp) and the second column %is the standard deviation. If the uncertainty is unknown, make the second %column all zero. % % INPUT tosave is the full path and name where the results will be % saved % TW = gapCloseParam.timeWindow % minTL = gapCloseParam.minTrackLen % LM = parameters.linearMotion % minSR = parameters.minSearchRadius % maxSR = parameters.maxSearchRadius 112 function error=testdetect(density,r, dm) % Generates a 3D image with spherical cells, samples it with 17 2D planes % and then tries to compute the centers of the cells with SpotDetector. % INPUT density is the number of cells % r is the radius of a cell (in pixels); assuming they’re all spherical % and have the same size % dm choose the detection method % - 0: my method, bottom up % - 1: my method, top down % - 2: adapted 3D method, with maximum % - 3: adapted 3D method, with center of max % - 4: detection with u-track % OUTPUT error(1) is the sum of the distance between each original cell B.4. testdetect B.4 testdetect % % % and the closest computed (by mySpotDetector) cell error(2) is the difference between the new density of cells (computed by mySpotDetector) and the original density %% Parameters bside=254; % number of pixels on each side of the 2D slices uheight=4; % distance (in pixels) between each pair of z slices nslices=17; % number of z slices height=(nslices-1)*uheight+2*r; imagen=uint8(zeros(bside,bside,nslices)); 113 %% Take 2D images for i=1:density % sweep all generated cells cx=C(i,1); cy=C(i,2); cz=C(i,3); % figure out which integer z are comprised in the sphere zstrt=max([ceil(cz-r),r]); zend=min([floor(cz+r),height-r]); % for each plane with integer z comprised in the sphere B.4. testdetect %% Generate cells (spheres) % rng(’shuffle’); %doesn’t work in MatlabR2010a C=rand(density,3); C(:,1:2)=0.499+bside*C(:,1:2); C(:,3)=height*C(:,3); end %rescale centers C(:,3)=(C(:,3)-r)/4+1; % add noise to original image % noise=uint8(64*randn(size(imagen))); % imagen=imagen+noise; B.4. testdetect for z=zstrt:zend if mod(z-r,4)==0 % check which integer z correspond to z-slice (ie were imaged) d=z-cz; % distance between center of sphere and current plane R=sqrt(r^2-d^2); % radius of the intersection of the sphere with the plane k=(z-r)/4+1; % corresponding number of z-slice % figure out which integer y are comprised in the circular intersection y0=max([ceil(cy-R),1]); y1=min([floor(cy+R),bside]); % for each integer y comprised in the intersection for y=y0:y1 sec=sqrt(R^2-(y-cy)^2); % x distance from y axis (of the circle) to circunference at height y x0=max([round(cx-sec),1]); x1=min([round(cx+sec),bside]); imagen(y,x0:x1,k)=128; end end end 114 %% Spot Detector -- Image denoising % for i=1:nslices % detectionMask = denoise(imagen(:,:,i)); % end 115 %% Sort results [~,ord2]=sort(C(:,3));C=C(ord2,:); B.4. testdetect %% Spot Detector 3D -- My Algorithm switch dm case 0 movieInfo=mySpotDetector(imagen,0); case 1 movieInfo=mySpotDetector(imagen,1); case 2 movieInfo=spotDetector3D(imagen,0); case 3 movieInfo=spotDetector3D(imagen,1); case 4 movieInfo=detecttr(imagen); end newdens=size(movieInfo.xCoord,1); output=zeros(newdens,3); output(:,1)=movieInfo.xCoord(:,1); output(:,2)=movieInfo.yCoord(:,1); output(:,3)=movieInfo.zCoord(:,1); [~,ord1]=sort(output(:,3));output=output(ord1,:); cd /home/monica/Thesisv2/simulations % this function needs to be improved in the % - find a better way to quantify error % - check the scaling when measuring error, % sense to rescale movieInfo instead of the % careful with the offsets (ie pixels start following ways: I think it would make more centers, but one needs to be at 0.5, not 0) B.5 testtrack 116 function tf=testtrack(no_of_tracks,no_of_steps,diffcoeff) % This function still needs a way of measuring the errors. % This function tests u-track on no_of_tracks simulated tracks which follow % Brownian motion with diffusion coefficient diffcoef. Tracks are followed B.5. testtrack %% Compute error dist=zeros(newdens,density); for i=1:newdens for j=1:density dist(i,j)=norm(output(i,:)-C(j,:)); end end [minbycol,midx]=min(dist,[],1); error(1)=sum(minbycol)/density; error(2)=newdens-density; % for no_of_steps seconds. curr=pwd; [prnt]=fileparts(curr); movieInfo=repmat(struct(’xCoord’,[],’yCoord’,[],’zCoord’,[],’amp’,[]),no_of_steps,1); %% Code for generating single-state random walks %% (i.e. pure diffusion with a single diffusion coefficient) displacements = normrnd(0, sqrt(diffcoeff*tau), [no_of_tracks,no_of_steps-1, dim]); %% Calculate position at each step startingpoint=cat(3,x*rand(no_of_tracks,1,2),z*rand(no_of_tracks,1,1)); position = cumsum(cat(2,startingpoint,displacements),2); 117 % This part ’registers’ only those positions of the simulated tracks which % are inside a box of size x by x by z (which represents the imaged area). % for i=1:no_of_steps % xidx=intersect(find(position(:,i,1)>=0),find(position(:,i,1)<=x)); % yidx=intersect(find(position(:,i,2)>=0),find(position(:,i,2)<=x)); % zidx=intersect(find(position(:,i,3)>=0),find(position(:,i,3)<=z)); B.5. testtrack %% Specify simulation parameters dim = 3; % Dimensionality tau = 6; % Sampling interval in s (0.001 = 1 ms) x= 250; z= 24; % % % % % % % end idx=intersect(xidx,yidx); idx=intersect(idx,zidx); movieInfo(i).xCoord=[position(idx,i,1),zeros(length(idx),1)]; movieInfo(i).yCoord=[position(idx,i,2),zeros(length(idx),1)]; movieInfo(i).zCoord=[position(idx,i,3),zeros(length(idx),1)]; movieInfo(i).amp=[ones(length(idx),1),zeros(length(idx),1)]; cd([prnt,filesep,’u-track110523’,filesep,’u-track’]) tf=scriptTrackGeneral(movieInfo,[prnt,filesep,’Results’,filesep,’testtr-’,... num2str(no_of_tracks),’-’,num2str(diffcoeff,’%1.2f’),’.mat’],dim); cd(prnt) B.5. testtrack % This part keeps full tracks even if they wander off the ’imaged’ box for i=1:no_of_steps movieInfo(i).xCoord=[position(:,i,1),zeros(no_of_tracks,1)]; movieInfo(i).yCoord=[position(:,i,2),zeros(no_of_tracks,1)]; movieInfo(i).zCoord=[position(:,i,3),zeros(no_of_tracks,1)]; movieInfo(i).amp=[ones(no_of_tracks,1),zeros(no_of_tracks,1)]; end 118 Appendix C Additional plots LN2sc − Step size LN2sc − Turning angle 16000 14000 5000 Frequency Frequency 12000 4000 3000 2000 10000 8000 6000 4000 1000 2000 0 0 Pi/4 Pi/2 3Pi/4 0 0 Pi 5 3500 300 3000 250 2500 Frequency Frequency 350 200 150 100 15 20 25 2000 1500 1000 50 0 0 10 Distance (um) LN2sc − Diameter Angle LN2sc − Average step size 500 2 4 6 8 10 0 0 12 50 Distance (um) 100 150 Distance (um) LN2sc − Track length 2000 Frequency 1500 1000 500 0 0 10 20 30 40 50 60 Time steps Figure C.1: Histograms for data set LN2 using mySpotDetector-bup as detection method with the lower level of denoising. 119 Appendix C. Additional plots LN8 − Turning angle LN8 − Step size 14000 12000 Frequency Frequency 2000 1500 1000 10000 8000 6000 4000 500 2000 0 0 Pi/4 Pi/2 3Pi/4 0 0 Pi 5 350 3000 300 2500 250 200 150 100 15 20 2000 1500 1000 500 50 0 0 10 Distance (pixels) LN8 − Diameter Frequency Frequency Angle LN8 − Average step size 2 4 6 8 0 0 10 20 Distance (pixels) 40 60 80 Distance (pixels) LN8 − Track length 3000 Frequency 2500 2000 1500 1000 500 0 0 10 20 30 40 50 60 Time steps Figure C.2: Histograms for data set LN8 using mySpotDetector-bup as detection method with the lower level of denoising. 120 Appendix C. Additional plots LN9 − Turning angle LN9 − Step size 7000 1200 6000 Frequency Frequency 1000 800 600 400 200 5000 4000 3000 2000 1000 0 Pi/4 Pi/2 3Pi/4 0 0 Pi 5 Angle LN9 − Average step size 10 15 20 Distance (pixels) LN9 − Diameter 300 2000 250 Frequency Frequency 1500 200 150 100 1000 500 50 0 0 2 4 6 8 10 0 0 12 10 Distance (pixels) 20 30 40 50 60 Distance (pixels) LN9 − Track length 1400 Frequency 1200 1000 800 600 400 200 0 0 5 10 15 20 Time steps Figure C.3: Histograms for data set LN9 using mySpotDetector-bup as detection method with the lower level of denoising. 121 Appendix C. Additional plots z1sc − Turning angle z1sc − Step size 5000 4000 Frequency Frequency 800 600 400 200 3000 2000 1000 0 Pi/4 Pi/2 3Pi/4 0 0 Pi 5 250 1000 200 800 150 100 50 0 0 10 15 20 Distance (um) z1sc − Diameter Frequency Frequency Angle z1sc − Average step size 600 400 200 2 4 6 8 0 0 10 10 Distance (um) 20 30 40 50 Distance (um) z1sc − Track length 2500 Frequency 2000 1500 1000 500 0 0 10 20 30 40 Time steps Figure C.4: Histograms for data set z1 using mySpotDetector-bup as detection method with the lower level of denoising. 122 Appendix C. Additional plots LN10 detection with utrack − angles LN10 detection with utrack − step size 200 25 150 Frequency Frequency 20 15 10 100 50 5 0 Pi/4 Pi/2 0 0 3Pi/4 2 Angle LN10 detection with utrack − diameter 6 8 10 LN10 detection with utrack − length 100 120 100 Frequency 80 Frequency 4 Distance (um) 60 40 80 60 40 20 0 0 20 5 10 15 Distance (um) 20 25 0 0 5 10 15 20 25 30 Time steps Figure C.5: Histograms corresponding to data set LN10 using u-track as method of detection, the z-coordinate determined as the normal average of all planes which contain the cell. For the tracking step, a pure Brownian motion model was assumed. 123 Appendix C. Additional plots LN10 detection with utrack − angles LN10 detection with utrack − step size 30 200 25 Frequency Frequency 150 20 15 100 10 50 5 0 Pi/4 Pi/2 0 0 3Pi/4 2 Angle 6 8 10 LN10 detection with utrack − length 140 100 120 100 80 Frequency Frequency LN10 detection with utrack − diameter 120 60 40 80 60 40 20 0 0 4 Distance (um) 20 5 10 15 Distance (um) 20 25 0 0 5 10 15 20 25 30 Time steps Figure C.6: Histograms corresponding to data set LN10 using u-track as method of detection, the z-coordinate determined as the normal average of all planes which contain the cell. For the tracking step, a random motion plus movement with constant velocity model was assumed. 124 Appendix C. Additional plots LN10 detection with utrack − angles LN10 detection with utrack − step size 30 250 200 20 Frequency Frequency 25 15 150 100 10 50 5 0 Pi/4 Pi/2 0 0 3Pi/4 2 4 Angle 8 10 12 LN10 detection with utrack − length 140 100 120 100 80 Frequency Frequency LN10 detection with utrack − diameter 120 60 40 80 60 40 20 0 0 6 Distance (um) 20 5 10 15 Distance (um) 20 25 0 0 5 10 15 20 25 30 Time steps Figure C.7: Histograms corresponding to data set LN10 using u-track as method of detection, the z-coordinate determined as the normal average of all planes which contain the cell. For the tracking step the motion model was random motion and movement along a straight line but with the possibility of immediate direction reversal. 125 Appendix C. Additional plots LN2trd − Step size LN2trd − Turning angle 4000 300 3500 250 Frequency Frequency 3000 200 150 100 2500 2000 1500 1000 50 0 500 Pi/4 Pi/2 0 0 3Pi/4 5 LN2trd − Diameter 15 20 LN2trd − Track length 2500 1200 1000 2000 800 Frequency Frequency 10 Distance (um) Angle 600 1500 1000 400 500 200 0 0 10 20 30 Distance (um) 40 50 0 0 20 40 60 80 100 120 Time steps Figure C.8: These histograms correspond to the analysis of the endogenous data set LN2, using u-track as method of detection and with the higher level of denoising. 126 Appendix C. Additional plots LN8trd − Turning angle LN8trd − Step size 700 50 600 500 Frequency Frequency 40 30 20 400 300 200 10 0 100 Pi/4 Pi/2 0 0 3Pi/4 5 Angle 10 15 20 Distance (um) LN8trd − Diameter LN8trd − Track length 250 400 350 300 Frequency Frequency 200 150 100 250 200 150 100 50 50 0 0 5 10 15 Distance (um) 20 25 0 0 10 20 30 40 50 60 Time steps Figure C.9: These histograms correspond to the analysis of the endogenous data set LN8, using u-track as method of detection and with the higher level of denoising. 127 Appendix C. Additional plots LN9trd − Turning angle LN9trd − Step size 200 150 Frequency Frequency 15 10 5 0 100 50 Pi/4 Pi/2 0 0 3Pi/4 2 4 Angle 70 140 60 120 50 40 30 80 60 40 10 20 10 Distance (um) 10 100 20 5 8 LN9trd − Track length 160 Frequency Frequency LN9trd − Diameter 80 0 0 6 Distance (um) 15 20 0 0 5 10 15 20 Time steps Figure C.10: These histograms correspond to the analysis of the endogenous data set LN9, using u-track as method of detection and with the higher level of denoising. 128 Appendix C. Additional plots z1trd − Step size 50 4 40 Frequency Frequency z1trd − Turning angle 5 3 2 1 0 30 20 10 Pi/4 Pi/2 0 0 3Pi/4 10 Angle z1trd − Diameter 40 50 40 Frequency 15 Frequency 30 z1trd − Track length 20 10 5 0 0 20 Distance (um) 30 20 10 10 20 Distance (um) 30 40 0 0 5 10 15 20 25 Time steps Figure C.11: These histograms correspond to the analysis of the endogenous data set Z1, using u-track as method of detection and with the higher level of denoising. 129
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Two mathematical approaches to a study of T cell motion...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Two mathematical approaches to a study of T cell motion and activation in the lymph node Delgado Carrillo, Monica 2012
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Two mathematical approaches to a study of T cell motion and activation in the lymph node |
Creator |
Delgado Carrillo, Monica |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | T cells are part of the immune system and as such play a very important role in keeping us healthy. One crucial step in the complex process which is the immune response to pathogens is T cell activation. The general goal of my thesis is to mathematically describe the migration patterns followed by T cells while waiting to be activated in the lymph node. Insight into these migration patterns could lead to better knowledge of the strategies T cells take to make activation such an efficient process. In order to fulfill my goal I have used two different approaches: one mainly computational and the other mainly theoretical. On the computational side, I analyzed three-dimensional microscopic movies of mice lymph nodes inside of which labelled T cells are moving. From the movies I extracted the trajectories of the cells. I studied movies from two experimental frameworks, exogenous and endogenous. On the former, more frequent type of experiment, T cells are labelled outside the mouse and then transferred in. The endogenous experiments, on the contrary, involve genetically modified mice whose T cells are born labelled. I concluded that there is a significant difference in labelled T cell motion between the two experimental frameworks. This suggests that previous results from exogenous experiments should be treated with caution due to possible errors introduced by the methods specific to that type of experiment. On the theoretical side I studied the time it takes for a model T cell to be activated under different scenarios regarding the characteristics of the lymph node as well as of the other cells in it. Since T cells become activated after establishing contact with a specific cell among many similar ones which also move within the lymph node, what I effectively computed was the mean first passage time for a model T cell to reach a defined target within the model lymph node. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-08-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial 3.0 Unported |
DOI | 10.14288/1.0073088 |
URI | http://hdl.handle.net/2429/43086 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2012-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc/3.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2012_fall_delgadocarrillo_monica.pdf [ 1.71MB ]
- Metadata
- JSON: 24-1.0073088.json
- JSON-LD: 24-1.0073088-ld.json
- RDF/XML (Pretty): 24-1.0073088-rdf.xml
- RDF/JSON: 24-1.0073088-rdf.json
- Turtle: 24-1.0073088-turtle.txt
- N-Triples: 24-1.0073088-rdf-ntriples.txt
- Original Record: 24-1.0073088-source.json
- Full Text
- 24-1.0073088-fulltext.txt
- Citation
- 24-1.0073088.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}]}"
data-media="{[{embed.selectedMedia}]}"
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-0073088/manifest