Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Two mathematical approaches to a study of T cell motion and activation in the lymph node Delgado Carrillo, Monica 2012

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
ubc_2012_fall_delgadocarrillo_monica.pdf [ 1.71MB ]
[if-you-see-this-DO-NOT-CLICK]
Metadata
JSON: 1.0073088.json
JSON-LD: 1.0073088+ld.json
RDF/XML (Pretty): 1.0073088.xml
RDF/JSON: 1.0073088+rdf.json
Turtle: 1.0073088+rdf-turtle.txt
N-Triples: 1.0073088+rdf-ntriples.txt
Original Record: 1.0073088 +original-record.json
Full Text
1.0073088.txt
Citation
1.0073088.ris

Full Text

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 FULFILLMENTOF 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 2012Abstract 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 migra- tion 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 con- trary, 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 re- sults from exogenous experiments should be treated with caution due to possible errors introduced by the methods specific to that type of experi- ment. 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. iiPreface My supervisor, Daniel Coombs has been involved in all areas of this the- sis. 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 Inflam- mation at the University of Glasgow, in Scotland, UK. Dr. Maffia and his collaborators own the copy rights to the data. All experiments were per- formed by Ross McQueenie. McQueenie also helped in describing the ex- perimental 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 Publish- ing 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. iiiTable of contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Biological description of the problem . . . . . . . . . . . . . 1 1.2 Current approaches . . . . . . . . . . . . . . . . . . . . . . . 2 2 Analysis of imaging data from T cells . . . . . . . . . . . . . . . 4 2.1 Two-photon microscopy . . . . . . . . . . . . . . . . . . . . . 4 2.2 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2.1 Experimental methods . . . . . . . . . . . . . . . . . 7 2.3 Methods for the inference of cell paths . . . . . . . . . . . . . 7 2.3.1 Single particle tracking . . . . . . . . . . . . . . . . . 7 2.3.2 Steps in the tracking process . . . . . . . . . . . . . . 9 2.3.3 Denoising with wavelets . . . . . . . . . . . . . . . . 10 2.3.4 Track extraction . . . . . . . . . . . . . . . . . . . . . 10 2.4 The detection step . . . . . . . . . . . . . . . . . . . . . . . . 14 2.4.1 Spot detector . . . . . . . . . . . . . . . . . . . . . . . 14 2.4.2 My detection technique . . . . . . . . . . . . . . . . . 15 2.4.3 Using u-track for detection . . . . . . . . . . . . . . . 17 ivTable of contents 2.5 Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.6 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.6.1 Discrimination of tracks . . . . . . . . . . . . . . . . . 25 2.6.2 Exogenous vs endogenous . . . . . . . . . . . . . . . 31 2.7 Discussion and future work . . . . . . . . . . . . . . . . . . . 34 3 First passage time of a T cell to an APC . . . . . . . . . . . . . . . 38 3.1 First passage time concepts . . . . . . . . . . . . . . . . . . . 39 3.2 Basic problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2.1 Neumann boundary . . . . . . . . . . . . . . . . . . . 40 3.2.2 Dirichlet boundary . . . . . . . . . . . . . . . . . . . . 43 3.2.3 Robin boundary . . . . . . . . . . . . . . . . . . . . . 49 3.3 Multiple targets with asymptotically null volume . . . . . . 54 3.3.1 Neumann boundary . . . . . . . . . . . . . . . . . . . 55 3.3.2 Robin boundary . . . . . . . . . . . . . . . . . . . . . 66 3.4 Space dependant diffusivity . . . . . . . . . . . . . . . . . . . 71 3.5 Discussion and future work . . . . . . . . . . . . . . . . . . . 75 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Appendices A Additional derivations . . . . . . . . . . . . . . . . . . . . . . . . 81 A.1 From the Neumann splitting probability, equation 3.51 . . . 81 A.1.1 The Neumann Green’s function for the sphere . . . . 81 A.1.2 Derivation of χ1 . . . . . . . . . . . . . . . . . . . . . 86 A.2 From the Robin splitting probability, equation 3.70 . . . . . 88 A.2.1 Derivation of χR . . . . . . . . . . . . . . . . . . . . . 88 B Programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 B.1 spotDetector3D . . . . . . . . . . . . . . . . . . . . . . . . . . 91 B.2 mySpotDetector . . . . . . . . . . . . . . . . . . . . . . . . . 99 B.3 detecttr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 B.4 testdetect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 B.5 testtrack . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 C Additional plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 vList of tables 2.1 Summary of the specifications of the experiments I analyzed data from . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Summary of the detection methods used . . . . . . . . . . . . 18 2.3 Number of tracks which are effectively moving in each ex- periment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.4 Maximum andminimum values of bootstrappedaverage ve- locities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.5 Confidence intervals of 95% of vexo − vendo . . . . . . . . . . . 34 viList of figures 1.1 Schematic diagram showing themajor structural components of a lymph node . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.1 Two photon microscopy requires two simultaneous events of high wavelength . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Two photon microscopy presents less background noise . . . 6 2.3 A sample microscopy image from each data set . . . . . . . . 8 2.4 A sample image in original and after denoising in the two levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.5 Sample trajectory computed by u-track . . . . . . . . . . . . 13 2.6 Extract from one of themovies to showhowmySpotDetector works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.7 Resulting saved regions after processing the example shown in Figure 2.6 with mySpotDetector . . . . . . . . . . . . . . . 17 2.8 Two 3D views of the connected region from the example in Figure 2.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.9 Sample simulated images . . . . . . . . . . . . . . . . . . . . 19 2.10 Error 1 as a function of the number of simulated cells . . . . 20 2.11 Error 2 as a function of the number of simulated cells . . . . 21 2.12 Histograms for an endogenousdata set usingmySpotDetector- bup as detection method with the lower level of denoising . 23 2.13 Histograms for an endogenousdata set using u-track as method of detection with the lower level of denoising . . . . . . . . . 24 2.14 Histograms for an endogenousdata set using u-track as method of detection with the higher level of denoising . . . . . . . . 25 2.15 Histograms for an exogenous using u-track as method of de- tection with the higher level of denoising . . . . . . . . . . . 26 2.16 Histogramof turning angle ignoring the z-coordinate for cho- sen tracks from all experiments with the higher level of de- noising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.17 3D representation of all tracks found for a data set . . . . . . 28 viiList of figures 2.18 Planar projections of Figure 2.17 . . . . . . . . . . . . . . . . 28 2.19 Turning angle histograms for odd and even time points sep- aratedly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.20 Example of how the process often misses part of the trajectory 30 2.21 Histograms of step size corresponding to the endogenous experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.22 Histograms of step size corresponding to the exogenous ex- periments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.23 Distributions of bootstrapped average cell velocity . . . . . . 33 2.24 Cells change shape when they are on motion . . . . . . . . . 35 2.25 Microscopy images showing the contour of the LN . . . . . . 36 2.26 A trajectory found by eye which the computational process missed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.1 Plots of the mean and variance of the FPT to the only target inner sphere with Neumann boundary condition . . . . . . . 42 3.2 Plot of the splitting probability to the only target inner sphere with Dirichlet boundary condition . . . . . . . . . . . . . . . 45 3.3 Plots of the mean and variance of the FPT to the only target inner sphere with Dirichlet boundary condition . . . . . . . 50 3.4 Plot of the splitting probability to the only target inner sphere with Robin boundary condition . . . . . . . . . . . . . . . . . 53 3.5 Plot of the splitting probability to the only target inner sphere with Robin boundary condition as a function of λ . . . . . . 53 3.6 Plots of the mean and variance of the FPT to the only target inner sphere with Robin boundary condition . . . . . . . . . 54 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 λ 55 3.8 Model of the LN with 3000 APCs . . . . . . . . . . . . . . . . 60 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 . . . . . . . . . . . . . . . . . . . . . . . 61 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 . . . . . . . . . . . . . . . . . . . . . . . 61 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 viiiList 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 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 . . . . . . . . . . 63 3.14 Plots of the space-averaged MFPT to the one specific target as a function of e and the number of targets when all targets have Dirichlet BC and the LN has Neumann BC . . . . . . . 63 3.15 Plots of the MFPT to the one specific target along two radii when all targets and the LN have Neumann BC . . . . . . . 66 3.16 Plots of the MFPT to the one specific target along two circles when all targets and the LN have Neumann BC . . . . . . . 67 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 . . . . . . . . . . . . . . . . . . . . . . . . . 70 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 . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.19 Two-dimensional version of the FRN model . . . . . . . . . . 71 3.20 A possible function to describe how the diffusion coefficient changes on the FRN . . . . . . . . . . . . . . . . . . . . . . . . 73 3.21 Comparison of theMFPT to the only target inner spherewith different BC on the LN . . . . . . . . . . . . . . . . . . . . . . 75 C.1 Histograms for data set LN2 using mySpotDetector-bup as detection method with the lower level of denoising . . . . . 119 C.2 Histograms for data set LN8 using mySpotDetector-bup as detection method with the lower level of denoising . . . . . 120 C.3 Histograms for data set LN9 using mySpotDetector-bup as detection method with the lower level of denoising . . . . . 121 C.4 Histograms for data set z1 using mySpotDetector-bup as de- tection method with the lower level of denoising . . . . . . . 122 C.5 Histograms corresponding to data set LN10 assuming a pure Brownian motion model for the tracking step . . . . . . . . . 123 C.6 Histograms corresponding to data set LN10 assuming a ran- dom motion plus movement with constant velocity model for the tracking step . . . . . . . . . . . . . . . . . . . . . . . . 124 ixList of figures C.7 Histograms corresponding to data set LN10 when the mo- tion model assumed for the tracking step was random mo- tion and movement along a straight line but with the possi- bility of immediate direction reversal . . . . . . . . . . . . . . 125 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 . . . . . . . . . . . . . . . . . . . 126 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 . . . . . . . . . . . . . . . . . . . 127 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 . . . . . . . . . . . . . . . . . . . 128 C.11 Histograms corresponding to the analysis of the endogenous data set Z1 using u-track as method of detection andwith the higher level of denoising . . . . . . . . . . . . . . . . . . . . . 129 xList 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 xiAcknowledgements I would like to thank my advisor, Daniel Coombs, who took me as his stu- dent 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 op- portunities 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 wel- coming me, and all my friends and family, here and in Mexico, who helped me through the good and bad moments. xiiDedication To my beloved Mexico. xiiiChapter 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 im- mune 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 patho- gen 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 im- portant 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 bloodstreamand migrate to the peripheral lymphoid organs (PLO), organized tissues where adaptive immune responses are initiated andwhere lymphocytes are main- tained. T cells are continually recirculating through these tissues. They en- ter 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 Pressent- ing 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 sys- 1Antigen - Term used to refer to any substance that can be recognized by the adaptive immune system. 11.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 net- work (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 equippedwith 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 prolifer- ation and differentiation, the first steps of adaptive immunity. The descen- dants 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 21.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 de- scribing 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 (two- photon microscopy) we were able to acquire time series of quasi-three- dimensional images showing T cells inside a LN. From these I extracted cell paths and characterized them. This has been done in the past, how- ever, 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. 3Chapter 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 im- age 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 specifica- tions are given in Table 2.1. 2.1 Two-photon microscopy Two-photon excitationmicroscopy [6] is a relatively new fluorescence imag- ing 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 en- ergy needed for single photon excitation, can equally excite a fluorophore if the two photons hit it simultaneously. Because that is a lower proba- bility event, fewer fluorophores will become excited, and more excitation events are likely to occur near the focal point, where more photons are con- centrated (Figure 2.1). This is in part the reason why 2PM presents less background noise (see Figure 2.2). Moreover, since lower energymeans longerwavelength, thewavelength of the absorbed photons in 2PM is about two times the wavelength the ab- sorbed 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] 42 .1 . T w o -ph oto n m icro scopy Name Exogenous/ Number of Number of Imaging rate Pixels µm Height per Endogenous slices time points (s/stack) per side per side slice (µm) LN2 endogenous 17 120 29.025 512 531.37 2 LN7 endogenous 12 150 20.004 512 531.37 2 LN8 endogenous 12 60 20.816 512 531.37 2 LN9 endogenous 30 20 50.799 512 531.37 2 LN10 exogenous 22 30 45.666 512 425.1 2 Z1 exogenous 25 38 12.821 560 283.4 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. 52.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 en- dogenous experiments on genetically modified mice whose cells express Green Fluorescent Protein (GFP). Because in these geneticallymodifiedmice 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 compli- 62.3. Methods for the inference of cell paths cated. 2.2.1 Experimental methods To image endogenous T cell populations, inguinal lymph nodes were ex- cised from hCD2 GFP mice and placed in a perfusion chamber. Lymph nodes were maintained at 37°C by perfusion with heated CO2 indepen- dent media. Following staining, approximately 2 million T cells were trans- ferred by intravenous tail vein injection into CD11c YFP transgenic mice. Multiphoton fluorescence was performed using a Zeiss 7MP multipho- ton microscope (Carl Zeiss, Jena, Germany) after excitation by a Coher- ent 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 sampleswere 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 mo- tion 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 dif- ficulties 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 par- ticle”. 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 72.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. 82.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 phenomenonwe 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 par- ticle 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. 92.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 Franc¸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 process- ing 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 res- olution 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 algo- rithm 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 oper- ations to get rid of isolated pixels and other outliers. This is done by means of the MATLAB function bwmorph, the higher level additionally perform- ing 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 denois- ing, unless otherwise stated. However, there are advantages and disad- vantages 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 102.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. 112.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, Brownianmotion plus movement with constant velocity, and randommotion plus movement along a straight line but with the pos- sibility of immediate direction reversal. I tested all three and the dif- ferences 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 us- ing the third motion model. It is worthmentioning 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 res- olution does not introduce bias to the results. To that effect the results of detection were transformed to micrometers using the resolution specifica- tions from each data set. 122.3. Methods for the inference of cell paths 152 154 156 158 116.5117117.5118118.5119119.5120120.5 15.5 16 16.5 17 17.5 18 18.5 19 19.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. 132.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 max- ima 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 re- turned 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 themax- imum 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. How- ever, I found that someone had already coded a function ordfilt3 [5] which nevertheless has the restriction of onlyworking 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 spotDetec- tor with ordfilt3 using a size 7 cubic window to define locality. Since scales 142.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 re- gion. 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 infor- mation. The idea behind my algorithm (hereafter referred to as mySpotDe- tector) 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 algo- rithm is explained next. 1. Find all 3D connected regions. I chose to do this defining connected- ness in terms of a 6-neighbors matrix, but I also did trials using 18- and 26-neighbors. 2. Loop over all these connected regions. 3. Look at the first (bottom most) plane that comprises the current re- gion. 4. Find all 2D connected regions and save them. For these I used an 8- neighbors 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 re- gions, append it to the list of saved regions. 152.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 tech- nique 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 re- gions, replace each of those saved regions with the correspond- ing intersection. 7. Repeat steps 5 and 6 until reaching the last plane that comprises the current 3D connected region. 8. Aweighted 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 (topmost) 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 162.4. The detection step “bottom up”. It does not appear to make much difference which ap- proach is used. Figures 2.6 to 2.8 show an example of how this detectionmethodworks. 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 Fig- ure 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 programwould have difficulty processing the information and often would show errors in 172.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 nor- mal 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 posi- tions 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 Section where Specific it is described parameter sd3D-MAX B.1. spotDetector3D 2.4.1 choice = 0 sd3D-COM B.1. spotDetector3D 2.4.1 choice = 1 mySD-bup B.2. mySpotDetector 2.4.2 bt = 0 mySD-tdown B.2. mySpotDetector 2.4.2 bt = 1 u-track B.3. detecttr 2.4.3 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. 182.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 detectionmethods 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 cor- responding 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 meth- ods. See Figure 2.9 for sample simulated planes. The interested reader can find the code for this test on Appendix B.4. 50 100 150 200 250 50 100 150 200 250 (a) 10 50 100 150 200 250 50 100 150 200 250 (b) 100 50 100 150 200 250 50 100 150 200 250 (c) 500 50 100 150 200 250 50 100 150 200 250 (d) 1000 50 100 150 200 250 50 100 150 200 250 (e) 2500 50 100 150 200 250 50 100 150 200 250 (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. 192.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, themethod that does best ismySpot- Detector. 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 test- ing 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 202.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 toomuch 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 212.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 mea- suring 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 parti- cle. 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. Ef- fectively 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 thenmade histograms showing this information. For example, Figures 2.12 and 2.14 both show results from data set LN7 but the first one corre- sponds 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 222.6. Analysis 0 Pi/4 Pi/2 3Pi/4 Pi0 1000 2000 3000 4000 5000 6000 LN7sc − Turning angle Angle Frequenc y (a) 0 20 40 60 80 0 0.5 1 1.5 2 2.5 3 3.5 x 10 4 LN7sc − Step size Velocity (um per minute) Frequenc y (b) 0 10 20 30 40 0 100 200 300 400 500 600 700 LN7sc − Average step size um per minute Frequenc y (c) 0 20 40 60 80 100 120 0 500 1000 1500 2000 2500 3000 LN7sc − Diameter Distance (um) Frequenc y (d) 0 50 100 150 0 1000 2000 3000 4000 5000 LN7sc − Track length Time steps Frequenc y (e) Figure 2.12: Histograms for data set LN7 using mySpotDetector-bup as de- tection method with the lower level of denoising. 232.6. Analysis Pi/4 Pi/2 3Pi/40 20 40 60 80 100 LN7trd − Turning angle Angle Frequenc y (a) 0 5 10 15 20 0 200 400 600 800 1000 1200 LN7trd − Step size Distance (um) Frequenc y (b) 0 5 10 15 20 0 50 100 150 200 250 LN7trd − Diameter Distance (um) Frequenc y (c) 0 20 40 60 80 100 0 200 400 600 800 1000 LN7trd − Track length Time steps Frequenc y (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. 242.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 pi would not be expected for random motion. In fact, no preferred angle would be expected for random motion. Pi/4 Pi/2 3Pi/40 10 20 30 40 50 60 70 80 LN7trd − Turning angle Angle Frequenc y (a) 0 5 10 15 20 0 100 200 300 400 500 600 700 800 LN7trd − Step size Distance (um) Frequenc y (b) 0 5 10 15 20 0 50 100 150 200 LN7trd − Diameter Distance (um) Frequenc y (c) 0 20 40 60 80 100 0 100 200 300 400 500 600 700 LN7trd − Track length Time steps Frequenc y (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 252.6. Analysis Pi/4 Pi/2 3Pi/40 5 10 15 20 25 LN10trd − Turning angle Angle Frequenc y (a) 0 2 4 6 8 10 12 0 50 100 150 200 250 300 LN10trd − Step size Distance (um) Frequenc y (b) 0 2 4 6 8 10 0 10 20 30 40 50 60 70 80 LN10trd − Average step size Distance (um) Frequenc y (c) 0 10 20 30 40 0 20 40 60 80 100 LN10trd − Diameter Distance (um) Frequenc y (d) 0 5 10 15 20 25 30 0 20 40 60 80 100 120 140 160 LN10trd − Track length Time steps Frequenc y (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. 262.6. Analysis 0 Pi/4 Pi/2 3Pi/4 Pi0 10 20 30 40 50 LN2trd − Turning angle ignoring z−coordinate Angle Frequenc y 0 Pi/4 Pi/2 3Pi/4 Pi0 10 20 30 40 50 LN7trd − Turning angle ignoring z−coordinate Angle Frequenc y 0 Pi/4 Pi/2 3Pi/4 Pi0 2 4 6 8 10 LN8trd − Turning angle ignoring z−coordinate Angle Frequenc y 0 Pi/4 Pi/2 3Pi/4 Pi0 0.5 1 1.5 2 LN9trd − Turning angle ignoring z−coordinate Angle Frequenc y 0 Pi/4 Pi/2 3Pi/4 Pi0 5 10 15 20 25 LN10trd − Turning angle ignoring z−coordinate Angle Frequenc y 0 Pi/4 Pi/2 3Pi/4 Pi0 1 2 3 4 5 6 Z1trd − Turning angle ignoring z−coordinate Angle Frequenc y 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). 272.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). 100 200 300 400 500 600 0 200 400 600 0 5 10 15 20 25 30 35 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. 100 150 200 250 300 350 400 450 500 550 0 100 200 300 400 500 600 (a) x - y 100 150 200 250 300 350 400 450 500 550 0 5 10 15 20 25 30 35 (b) x - z 0 100 200 300 400 500 600 0 5 10 15 20 25 30 35 (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 histogramswould show the true turning angle distribution. This 282.6. Analysis Pi/4 Pi/2 3Pi/40 5 10 15 20 25 30 35 LN8trd − Turning angle skipping even positions Angle Frequenc y Pi/4 Pi/2 3Pi/40 5 10 15 20 25 LN8trd − Turning angle skipping odd positions Angle Frequenc y Figure 2.19: Turning angle histograms for odd and even time points sepa- ratedly. 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 observedmany 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 trav- elled during a time step is less than 1µm for tracks that correspond to sta- tionary cells. For that reason I computed the number of tracks whose av- erage 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. 292.6. Analysis Name Level of denoising Tracks w/ big avg steps Long tracks LN2 higher 540 100 LN2 normal 753 140 LN7 higher 326 56 LN7 normal 530 90 LN8 higher 161 30 LN8 normal 248 34 LN9 higher 23 3 LN9 normal 43 9 LN10 higher 98 40 LN10 normal 111 41 Z1 higher 76 14 Z1 normal 158 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 tra- jectory. 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. 302.6. Analysis 0 10 20 30 40 50 60 0 1000 2000 3000 4000 5000 6000 LN2trd − Step size Velocity (um per minute) Frequenc y 0 10 20 30 40 50 60 0 200 400 600 800 1000 LN7trd − Step size Velocity (um per minute) Frequenc y 0 10 20 30 40 50 60 0 100 200 300 400 500 600 700 LN8trd − Step size Velocity (um per minute) Frequenc y 0 10 20 30 40 50 60 0 100 200 300 400 500 600 LN9trd − Step size Velocity (um per minute) Frequenc y Figure 2.21: Histograms of step size corresponding to the endogenous ex- periments. 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 312.6. Analysis 0 10 20 30 40 50 60 0 50 100 150 200 250 300 350 LN10trd − Step size Velocity (um per minute) Frequenc y 0 50 100 150 200 0 5 10 15 20 Z1trd − Step size Velocity (um per minute) Frequenc y Figure 2.22: Histograms of step size corresponding to the exogenous ex- periments. 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 endoge- nous experiments and 1,349 data points from the exogenous experiments. To obtain the distribution of mean velocity from this data, I took 1000 boot- strap 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 experi- ment 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 experimentsmove significantly 322.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 re- sult is highly suggestive of an important difference in labelled T cell motion between exogenous and endogenous experiments. Similar results were obtainedwhen using the normal level of denoising. Variable Minimum Maximum vendo 2.0902 2.2662 vexo 3.7452 5.0438 vZ1 11.0362 16.1753 vLN10 2.5323 2.9707 Table 2.4: Maximum and minimum values of bootstrapped average veloc- ities. Units are µm/min. vendo is the bootstrapped average of 1000 sam- ples collected from all four endogenous experiments while vexo is the cor- responding 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. 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Velocity (um per minute) Probabilit y   Figure 2.23: Distributions of bootstrapped average cell velocity: vendo in red, vexo in black, vZ1 in magenta, and vLN10 in blue. 332.7. Discussion and future work Variable Confidence interval vexo − vendo 1.6772 - 2.3442 vZ1 − vendo 9.3529 - 12.9622 vLN10 − vendo 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 be- cause 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 pi, would be overcome. There is another piece of information that I think should be incorpo- rated 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 342.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 imageswhere the bound- ary 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 bound- ary 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 sub- stracted 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 dis- placement showed. Finding a better way to measure the drift is on the list of future work. These and other factors are the reason whywe think the tracking results need a lot of improvement yet. Lack of adequate visualization tools is one more problem that, if over- come, 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 soft- ware 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 plan- ning 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 352.7. Discussion and future work Figure 2.25: Microscopy images showing the contour of the LN. These im- ages correspond to the first and last time points, plane z = 1, from data set LN7. 362.7. Discussion and future work Figure 2.26: This figure shows a trajectory found by eye which the compu- tational 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 vi- sualization of biological images. Unfortunately, I did not find a satisfactory solution. Besides improving the detection and tracking methods one natural ex- tension to this work is to analyze more experiments, especially of the ex- ogenous 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 move- ment patterns and zones in the LN, the influence of zones being not only suggested by theory but also by observing the videos. 37Chapter 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 big- ger 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 simu- lated 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 com- pared 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. Neverthe- less, cells do not exit when they first “bump into” this hole, but they con- 383.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 con- centration u (x, t) at time t on the point x in space, satisfies the diffusion equation ∂u ∂t = ∂ ∂x [ D (x) ∂u∂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 ∂nu (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 ∂nu (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 disjoint sets B = ∪Ni=1Bi, Pi is the probability of exiting through Bi. Moments of the first passage time Tn is the n-th moment of the first pas- sage 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. 393.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 al- ready 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 ∇2T (x) = − 1D , x ∈ ΩΩ1; (3.1)T (x) = 0, x ∈ ∂Ω1; ∂nT (x) = 0, x ∈ ∂Ω, 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 ∂2T ∂r2 (r) + 2 r ∂T ∂r (r) = − 1 D , r ∈ (R−, R+), (3.2)T (R − ) = 0, ∂nT (R+) = 0. The solution to the homogeneous version of the Poisson equation, that is the Laplace equation, is c1 r + c2, (3.3) 403.2. Basic problem and we find a particular solution for 3.2 to be − r26D . Hence the general solution to 3.2 is T (r) = c1 r + c2 − r2 6D . (3.4) Now the boundary conditions T (R − ) = 0 and ∂nT (R+) = 0 give the following equations for c1 and c2 c1 R − + c2 − R2 − 6D = 0, T′ (R+) = − c1R2+ − R+ 3D = 0, from where c1 =− R3+ 3D , c2 = R2 − 6D − c1 R − = R2 − 6D + R3+ 3DR − = R3 − + 2R3+ 6DR − . The final expression for T is T (r) = R 3 − + 2R3+ 6DR − − R3+ 3Dr − r2 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 ∇2T2 (x) = − 2DT (x) , x ∈ Ω; (3.7)T2 (x) = 0, x ∈ ΩΩ1; ∂nT2 (x) = 0, x ∈ ∂Ω. I found the following to be a particular solution for the previous equation: r4 60D2 + R3+r 3D2 − ( R3 − + 2R3+ ) r2 18D2R − . 413.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 T2 (r) = c1 r + c2 + r4 60D2 + R3+r 3D2 − ( R3 − + 2R3+ ) r2 18D2R − . After substituting the boundary conditions and solving for c1 and c2 we find a final expression for the second moment of the FPT (3.8) T2 (r) = 18R 5 +R− − 5R3+ ( R3 − + 2R3+ ) 45D2R − r + r ( 3R − r3 + 60R3+R− − 10r ( R3 − + 2R3+ )) 180D2R − + 40R6+ + 7R6− − 72R−R5+ − 20R3−R3+ 180D2R2 − , and so we also have an expression for the VFPT, which we denote by V (3.9) V (r) = r ( 20R3+ − r3 ) 90D2 + 10R6+ + R6− − 20R3−R3+ − 36R−R5+ 90D2R2 − + 2R5+ 5D2r − R6+ 9D2r2 . A plot ofV is shown on Figure 3.1. 423.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 bound- ary 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 prob- ability 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: ∇2P (x) = 0, x ∈ ΩΩ1, (3.10)P (x) = 0, x ∈ ∂Ω, P (x) = 1, x ∈ ∂Ω1; (3.11)D∇2 [PT] (x) = −P (x) , x ∈ ΩΩ1, [PT] (x) = 0, x ∈ ∂Ω ∪ ∂Ω1. Again by spherical symmetry the previous equations are reduced to P ′′ (r) + 2 r P ′ (r) = 0, (3.12)P (R − ) = 1, P (R+) = 0; D [ P ′′T + 2P ′T′ + PT′′ + 2 r ( P ′T + PT′ )] = −P , (3.13)P (R − )T (R − ) = 0, P (R+)T (R+) = 0. 433.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. P (R − ) = c1R − + c2 = 1 P (R+) = c1R+ + c2 = 0 =⇒ c2 = − c1 R+ =⇒ R+c1 − R−c1 = R−R+ =⇒ c1 = R − R+ R+ − R− =⇒ c2 = − R − R+ − R− Hence P (r) = R−R+ r (R+ − R−) − R − R+ − R− . (3.14) Equation 3.13 then reduces to T′′ + 2c2 c1 + c2r T′ = − 1D (3.15) and we have that 2c2 c1 + c2r = − 2 R+ − r , so the solution is T (r) = − r (r− 2R+)6D − k1 r− R+ + k2. Now we need to find k1 and k2 from the boundary conditions. We have P (r)T (r) = c2R+ (r − 2R+)6D − c2r (r − 2R+) 6D + k1c2R+ r (r − R+) − k1c2 r − R+ − k2c2R+ r + k2c2 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 1r−R+ term. We have 0 = P (R+)T (R+) = k1c2 (R+ − R+)R+ (R+ − R+) = − k1c2 R+ 443.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 − ) = −R− (R− − 2R+)6D + k2, and thus k2 = R − (R − − 2R+) 6D . The final expresion for T is T (r) = − r (r− 2R+)6D + R − (R − − 2R+) 6D . (3.16) Plots of P and T for the case with absorbing boundary are shown on Figures 3.2 and 3.3. Variance of the first passage time In [22], the author does not present equations for the higher ordermoments of the FPT. However, it is easy to deduce the equation for the second mo- ment 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. P (x) = ∑ p∈paths Pp (x) . (3.17) 453.2. Basic problem Then, by their definitions, T and T2 can be written as T (x) = ∑p∈paths Pp (x) tp (x)∑p∈paths Pp (x) , (3.18) T2 (x) = ∑p∈paths Pp (x) t 2 p (x) ∑p∈paths Pp (x) , (3.19) where tp (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 in- termediate 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 direc- tions (±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 P (x) =∑ p [ 3 ∑ i=1 1 6Pp (x + δxei) + 3 ∑ i=1 1 6Pp (x − δxei) ] = 1 6 3 ∑ i=1 [ Pp (x + δxei) + Pp (x − δxei) ] . 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 463.2. Basic problem P (x)T2 (x) = ∑ p∈paths Pp (x) t2p (x) =∑ p 3 ∑ i=1 [ 1 6Pp (x + δxei) ( tp (x + δxei) + δt )2 + 1 6Pp (x − δxei) ( tp (x − δxei) + δt )2] = δt2P (x) + δt3 3 ∑ i=1 [ ∑ p Pp (x + δxei) tp (x + δxei) +∑ p Pp (x − δxei) tp (x − δxei) ] + 1 6 3 ∑ i=1 [ ∑ p Pp (x + δxei) t2p (x + δxei) +∑ p Pp (x − δxei) t2p (x − δxei) ] = δt2P (x) + δt 3 3 ∑ i=1 [P (x + δxei)T (x + δxei) + P (x − δxei)T (x − δxei)] + 1 6 3 ∑ i=1 [P (x + δxei)T2 (x + δxei) + P (x − δxei)T2 (x − δxei)] where the last expression is obtained using equations 3.18 and 3.19. Next we replace [PTn] (x± δxei) by its three term Taylor expansion about x. Then P (x)T2 (x) = δt2P (x) + δt3 3 ∑ i=1 [ 2P (x)T (x) + δx2 ∂ 2PT ∂x2i (x) ] + 1 3 3 ∑ i=1 [ P (x)T2 (x) + δx 2 2 ∂2PT2 ∂x2i (x) ] = δt2P (x) + 2δtP (x)T (x) + δtδx 2 3 ∇ 2 [PT] (x) + P (x)T2 (x) + δx 2 2 ∇ 2 [PT2] (x) (3.20) =⇒ −2δtP (x)T (x) = δt2P (x) + δtδx 2 3 ∇ 2 [PT] (x) + δx 2 2 ∇ 2 [PT2] (x) . 473.2. Basic problem Now in this discretized diffusion process a relation between the time and space steps is enforced, namely that D = δx 2 2δt . Thus equation 3.20 turns into ( δx2 2D )2 P (x) + δx 4 6D∇ 2 [PT] (x) + δx 2 2 ∇ 2 [PT2] (x) = −δx 2 D P (x)T (x) . Finally, we divide the previous equation by δx2 and then take the limit as δx goes to 0 to obtain ∇2 [PT2] (x) = − 2DP (x)T (x) . (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 T′′2 + 2 r− R+ T′2 = 1 3D2 (r (r− 2R+) − R− (R− − 2R+)) (3.22) which has solution given by T2 (r) = 3D 2b1 + 5R3+R− (R− − 2R+) + 2R5+ 45D2 (r − R+) − (R − − R+)2 (r − R+)2 9D2 + (r − R+)4 20D2 + b2, where b1 and b2 are constants to be determined upon sustitution of bound- ary 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 1r−R − . Cancelling terms where applicable we are left with the following equation at r = R+ 0 = [PT2] (r) = R+ ( 3D2b1 + 5R3+R− (R− − 2R+) + 2R5+ ) − r ( 3D2b1 + 5R3+R− (R− − 2R+) + 2R5+ ) 45D2r (r − R+) = − (r − R+) 45D2r (r − R+) ( 3D2b1 + 5R3+R− (R− − 2R+) + 2R5+ ) . 483.2. Basic problem Hence b1 = − 5R3+R− (R− − 2R+) + 2R5+ 3D2 . Now, since P (R − ) = 1 at the other boundary we have 0 = T2 (R−) = − (R − − R+)4 9D2 + (R − − R+)4 20D2 + b2 and thus b2 = 11 (R − − R+)4 180D2 . Finally, we arrive to an expression for T2 T2 (r) = ( 11 (R − − R+)2 − 9 (r− R+)2 ) ( (R − − R+)2 − (r− R+)2 ) 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: (3.24) V (r) = ( 11 (R − − R+)2 − 9 (r − R+)2 ) ( (R − − R+)2 − (r − R+)2 ) 180D2 − ( (R − − R+)2 − (r − R+)2 ) 2 36D2 = ( (R − − R+)2 − (r − R+)2 ) ( 3 (R − − R+)2 − 2 (r − R+)2 ) 90D2 . 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 493.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 (3.25)P (R+) = (1− λ)P (R+ − δx) ≈ (1− λ) (P (R+)− δxP ′ (R+)) , 503.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+) = −λ0P (R+) . Taking the limit as δx goes to 0 we get the following condition for the split- ting probability on the right hand side boundary: P ′ (R+) = −λ0P (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 [PT] (R+) =∑ p (1− λ) Pp (R+ − δx) ( tp (R+ − δx) + δt ) = δt (1− λ)P (R+ − δx) + (1− λ)P (R+ − δx)T (R+ − δx) = δtP (R+) + (1− λ) ([PT] (R+)− δx [PT] ′ (R+)) , where we have subsituted equation 3.25 and again used a Taylor approxi- mation of PT about R+ to obtain the last line. This yields (1− λ) δx [PT]′ (R+) = δtP (R+) − λ [PT] (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 equa- tion and then dividing by δx we get (1− λ0δx) [PT]′ (R+) = δx2DP (R+) − λ0 [PT] (R+) . Again taking the limit as δx goes to 0 gives the boundary conditionwewere looking for: [PT]′ (R+) = −λ0 [PT] (R+) . (3.27) The same condition is satisfied by the second moment, as can be seen following the analogous argument, i.e. [PT2] (R+) =∑ p (1− λ) Pp (R+ − δx) ( tp (R+ − δx) + δt )2 = δt2P (R+) + 2δt (1− λ) [PT2] (R+ − δx) + (1− λ) ([PT2] (R+)− δx [PT2] ′ (R+)) from where [PT2]′ (R+) = −λ0 [PT2] (R+) . (3.28) 513.2. Basic problem Splitting probability and mean first passage time Nowwe 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 ab- sorbing 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 solu- tion to 3.10, namely P (r) = c1 r + 1− c1R − , Then condition 3.26 gives − c1 R2+ + λ c1R+ + λ ( 1− c1R − ) = 0, (3.29) which then implies −c1R− + λc1R+R− + λR2+R− − λc1R2+ = 0 and so c1 = λR2+R− λR2+ − λR+R− + R− , (3.30) c2 = R − (1− λR+) λR2+ − λR+R− + R− . (3.31) Hence the solution for P is P (r) = R− ( λR2+ − λR+r + 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′′ + 2 (1− λR+)λR2+ − λR+r + rT ′ = − 1 D . The general solution of the previous equation is T (r) = k2 − ( λR2+ − λR+r + r )2 6D (λR+ − 1)2 − k1 6D (λR+ − 1) (λR2+ − λR+r + r) . (3.33) 523.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+. 533.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 [PT]′ (R+) = −λ [PT] (R+), we obtain an expression for the MFPT: T (r) = 13D (λR+ − 1) [ ( λR2+ − λR+R− + R− )2 − ( λR2+ − λR+r + r )2 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 Ωei with centers xi ∈ Ω with radii eai, ai > 0, in such a way that these spheres are completely contained inside 543.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 rea- sons 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 Ωe1 before reaching any of the other APCs Ωej for j = 2, . . . , N, by using matched asymptotics. The starting point x is in Ω\Ωp where Ωp ≡ ∪Ni=1Ωei . I followed notes from 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 553.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 ∇2P (x) = 0, x ∈ Ω\Ωp; (3.35)∂nP (x) = 0, x ∈ ∂Ω; P (x) = 1, x ∈ ∂Ω1; P (x) = 0, x ∈ ∪Nj=2∂Ωei . We first look at the outer region, away from the targets. As usual, we expand P as P = P0 + eP1 + e2P2 + . . . (3.36) Then (3.37)∇2Pk = 0, x ∈ Ω\ {x1, . . . , xN} ∂nPk = 0, x ∈ ∂Ω with certain singularity conditions as x → xj 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 ( xj + ey ) , with y ≡ e−1 ( x− xj ) , w = w0 + ew1 + . . . (3.38) Now ∇2yw = e2∇2xP and so ∇2yw0 = 0, y /∈ Ωj; w0 = δj1, y ∈ ∂Ωj (3.39) ∇2yw1 = 0, y /∈ Ωj; w1 = 0, y ∈ ∂Ωj (3.40) Here Ωj = e−1Ωej . The far-field boundary conditions for w0 and w1 are determined by the matching condition as x → xj between the inner and outer expansions 3.36 and 3.38, written as P0 + eP1 + e2P2 + . . . ∼ w0 + ew1 + . . . (3.41) The first matching condition is that w0 ∼ P0 as |y| → ∞. As before, the solution to 3.39 is w0 = c1 |y| + c2, 563.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 ( aj ) = δj1, and so c1 = aj ( δj1 −P0 ) . Hence w0 (y) ≈ P0 + aj ( δj1 −P0 ) |y| (3.42) as |y| → ∞. From this and the matching condition 3.41 we conclude that as x → xj, P1 ∼ aj ( δj1 −P0 ) |x− xj| . (3.43) To be able to satisfy this condition we slightly perturb equation 3.37: ∇2P1 (x) = Sδ (x) , where S is a constant and δ is the Dirac delta function. Transforming the previous equation to spherical coordinates we have 1 r2 ∂ ∂r [ r2 ∂P1 ∂r ] = Sδ (r) , which upon integration holds ∂ ∂r [ r2 ∂P1 ∂r ] = Sr2δ (r) . (3.44) Now recall that by definition∫ R3 δ (x) dx = 1. We can transform this equation to spherical coordinates to obtain (3.45) 1 = ∫ 2pi 0 ∫ pi 0 ∫ ∞ 0 δ (r) r2 sin θdrdθdϕ = ∫ ∞ 0 δ (r) r2 ∫ 2pi 0 ∫ pi 0 sin θdθdϕdr = ∫ ∞ 0 δ (r) r2 ∫ 2pi 0 [− cos θ]pi0 dϕdr = ∫ ∞ 0 δ (r) r2 ∫ 2pi 0 2dϕdr = 4pi ∫ ∞ 0 δ (r) r2dr. 573.3. Multiple targets with asymptotically null volume Hence, using 3.45 on equation 3.44 we get r2 ∂P1 ∂r = S 4pi + c ′ 1 =⇒ ∂P1 ∂r = S 4pir2 + c′1 r2 =⇒ P1 (r) = −S + c ′′ 1 4pir + c ′ 2. Nowwe can find out what S has to be to agree with the matching condition 3.43, namely S = −4pi ( δj1 −P0 ) aj. Therefore the equation for P1 writes (3.46)∇2P1 = −4pi N ∑ j=1 ( δj1 − P0 ) ajδ ( x − xj ) , x ∈ Ω; ∂nP1 = 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∫∫∫ Ω ∇2P1dV = ∫∫ ∂Ω ∇P1 ·~ndS = ∫∫ ∂Ω ∂nP1dS and the RHS is equal to zero by the boundary condition on 3.46. On the other hand ∫∫∫ Ω ∇2P1dV = −4pi N ∑ j=1 ( δj1 −P0 ) aj and so (1−P0) a1 −P0∑Nj=2 aj = 0 (3.47) =⇒ P0 = a1 ∑Nj=1 aj . We write (3.48)P0 = a1Na¯ , a¯ ≡ a1 + · · · + aN N , to save space. 583.3. Multiple targets with asymptotically null volume Now let us look at the Green’s function for the system ∇2G = 1|Ω| − δ (x − ξ) , x ∈ Ω; (3.49)∂nG = 0, x ∈ ∂Ω;∫ Ω G (x; ξ) dx = 0. In terms of this function the solution to 3.46 can be written as P1 (x) = 4pi N ∑ j=1 ( δj1 −P0 ) ajG ( x; xj ) + χ1 (3.50) where χ1 is a constant. By integrating the previous equation we obtain an expression for χ1:∫ Ω P1dx = 4pi N ∑ j=1 ( δj1 − P0 ) aj ∫ Ω G ( x; xj ) + ∫ Ω χ1 = χ1|Ω| by equation 3.47. Then χ1 = 1 |Ω| ∫ Ω P1dx We have finally arrived to a two term expansion of the splitting proba- bility: P (x) = a1Na¯ + 4piea1 [ G (x; x1) − 1Na¯ N ∑ j=1 ajG ( x; xj )] + eχ1 +O ( e2 ) . (3.51) 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 ∇2TG (x) = − 1D , x ∈ Ω\Ωp, (3.52)∂nTG (x) = 0, x ∈ ∂Ω, TG (x) = 0, x ∈ ∂Ωp, 593.3. Multiple targets with asymptotically null volume −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Figure 3.8: Sphere of radius 1 representing the lymph node, with 3000 ran- domly 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 5e 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 aj = 0.35 for all j, since for e = 0.01 the radius eaj = 3.5× 10−3mm would correspond to the average size of a dendritic cell. However, such a big value for e made it imposible to compute posi- tions for the 3000 targets without overlapping so I used e = 0.001 instead. See the following figures for more details. 603.3. Multiple targets with asymptotically null volume 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 x 10 −3 r P 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 x 10 −3 r P Figure 3.9: Splitting probability from equations 3.51, A.11, and A.15, with e = 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). 0 1.5708 3.1416 4.7124 6.2832 1.5 2 2.5 3 3.5 4 4.5 5 5.5 x 10 −4 Azimuthal angle P 0 1.5708 3.1416 4.7124 6.2832 1 2 3 4 5 6 7 8 x 10 −4 Azimuthal angle P Figure 3.10: Splitting probability from equations 3.51, A.11, and A.15, with e = 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). 613.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. ThusTG is given in the limit e → 0 of small trap radius by (3.53)TG ≈ a 3 3Na¯De [ 1− 4pie N ∑ j=1 ajG ( x; xj ) + 4pie Na¯ a TGa +O (e2) ] [4]. See Figures 3.11 and 3.12 for examples. 0 0.2 0.4 0.6 0.8 1 0.086 0.088 0.09 0.092 0.094 0.096 0.098 0.1 r T 0 0.2 0.4 0.6 0.8 1 0.085 0.09 0.095 0.1 0.105 0.11 0.115 r T Figure 3.11: MFPT from equations 3.51 and A.11, with e = 0.001 and D = 3. TG 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 1.5708 3.1416 4.7124 6.2832 0.1 0.1005 0.101 0.1015 0.102 0.1025 0.103 0.1035 Azimuthal angle T 0 1.5708 3.1416 4.7124 6.2832 0.1 0.101 0.102 0.103 0.104 0.105 0.106 Azimuthal angle T Figure 3.12: MFPT from equations 3.51 and A.11, with e = 0.001 and D = 3. TG 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). 623.3. Multiple targets with asymptotically null volume Also, the average MFPT based on a uniform distribution of starting points, ¯TG = |Ω|−1 ∫ΩTGdx, satisfies (3.54)¯TG ≈ a 3 3Na¯De [ 1 + 4pieNa¯ a TGa +O (e2) ] . ¯TG for parameters like those from Figure 3.11 is around 0.1075. Increasing the radii of the targets to 3.5 gives a ¯TG of about 0.0125, while reducing e to 0.0001 makes ¯TG = 1.0885. See also Figures 3.13 and 3.14. 0.02 0.04 0.06 0.08 0.1 0 0.5 1 1.5 2 2.5 3 3.5 4 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 2 4 6 8 10 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 Figure 3.13: Plot of ¯TG versus aj. N = 3000, a = 1, e = 0.001 and D = 3. Due to the big variation in scale I separated the range of aj in three plots. 2 4 6 8 10 x 10−4 0.2 0.4 0.6 0.8 1 ε 500 1000 1500 2000 2500 3000 3500 4000 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 N Figure 3.14: Plot of ¯TG versus e 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 ∇2 [PT] (x) = − 1DP (x) , x ∈ Ω\Ωp; (3.55)∂nT (x) = 0, x ∈ ∂Ω; [PT] (x) = 0, x ∈ ∂ ∪Nj=2 Ωej ; T (x) = 0, x ∈ ∂Ωe1 . 633.3. Multiple targets with asymptotically null volume Unfortunately, we have not as yet been able to solve this. Instead, we study the previouslymentioned 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 infinites- imaly short. The differential equation to solve is ∇2T (x) = − 1D , x ∈ Ω\Ωp; (3.56)∂nT (x) = 0, x ∈ ∂Ω⋃∪Nj=2Ωej ; T (x) = 0, x ∈ ∂Ωe1 . In the outer region we expand T = e−1T0 +T1 + eT2 + . . .. For k = 0, 2 the Tk satisfy ∇2Tk (x) = 0, x ∈ Ω\ {x1, . . . , xN}, (3.57)∂Tk (x) = 0, x ∈ ∂Ω, while T1 satisfies ∇2T1 (x) = − 1D , x ∈ Ω\ {x1, . . . , xN}; (3.58)∂nT1 (x) = 0, x ∈ ∂Ω. Singularity conditions as x → xj will be determined upon matching to the inner solution. In the inner region near the j-th trap, like we did for P before, we let v (y) ≡ T (xj + ey) with y ≡ e−1 (x− xj). We also expand v = e−1v0 + v1 + . . .. The equations for the vk are ∇2vk ( y ) = 0, y /∈ Ωj, (3.59)∂nvk ( y ) = 0, y ∈ ∂Ωj, j 6= 1, vk ( y ) = 0, y ∈ ∂Ω1, and the matching condition reads e−1T0 +T1 + eT2 + . . . ∼ e−1v0 + 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) = T0 − T0a1|y| 643.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 ∼ − a1T0 |x− x1| as x → x1 and T1 ∼ 0 as x → xj. Therefore, in terms of the Dirac distribu- tion, T1 satisfies ∇2T1 (x) = 4pia1T0δ (x − x1)− 1D , x ∈ Ω; ∂nT1 (x) = 0, x ∈ ∂Ω. We use the divergence theorem on the previous equation to derive a solv- ability condition for T1, which yields T0 = |Ω| 4pia1D . (3.60) Now in terms of the Green’s function, defined by equation 3.49, and a constant χT we find that T1 = −4pia1T0G (x; x1) + χT . Next we expand the previous expression for T1 as x → xj using the local behaviour of G found on Appendix A.1.1. We get T1 ∼ − a1T0 |x − x1| − 4pia1T0R1,1 + χT, as x → x1; T1 ∼ −4pia1T0Gj,1 + χT, as x → xj. Upon substituting into the matching condition we find that v1 ∼ −4pia1T0Gj,1 + χT, as |y| → ∞, where we are using the matrix notation from Appendix A.1.1. From equa- tion 3.59 we note that v1 is also a constant in the inner region near the j-th trap for j 6= 1, whereas v1 (y) = c1 ( 1− a1|y| ) in the inner region near the first trap, with c1 a constant. Hence v1 (y) = (χT − 4pia1T0Gj,1) ( 1− δj,1a1 |y| ) . 653.3. Multiple targets with asymptotically null volume 0 0.2 0.4 0.6 0.8 1 317 317.05 317.1 317.15 317.2 317.25 317.3 317.35 317.4 r T 0 0.2 0.4 0.6 0.8 1 316.5 316.6 316.7 316.8 316.9 317 317.1 317.2 317.3 317.4 r T Figure 3.15: MFPT from equation 3.61 with e = 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 → xj, j 6= 1, and T2 ∼ a1 (4pia1T0R1,1 − χT) |x− x1| as x → x1. Thus the equation for T2 is ∇2T2 (x) = −4pi ( 4pia1T0R1,1 − χT ) a1δ (x − x1) , x ∈ Ω; ∂nT2 (x) = 0, x ∈ ∂Ω. Finally, the divergence theorem yields a solvability condition for T2 which gives χT = 4pia1T0R1,1. We have thus obtained a two-term expression for the MFPT: T (x) = a 3 3ea1D + 4pia3 3D (R1,1 − G (x; x1)) . (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 ∂Ω, ∂nP = 0, is replaced by ∂nP = λP . Since the Robin Green’s function for the Laplace equation is unknown, we instead consider the case when λ = eλ0 for some constant λ0 and let e 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. 663.3. Multiple targets with asymptotically null volume 0 1.5708 3.1416 4.7124 6.2832 317.2 317.22 317.24 317.26 317.28 317.3 317.32 317.34 Azimuthal angle T 0 1.5708 3.1416 4.7124 6.2832 317.12 317.14 317.16 317.18 317.2 317.22 317.24 317.26 317.28 317.3 Azimuthal angle T Figure 3.16: MFPT from equation 3.61 with e = 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 + eP1 + e2P2 + . . .. Notice that the boundary condition for P0 is again ∂nP0 = 0 so that P0 is a constant. Now, near the j-th trap we again let w (y) ≡ P (xj + ey) and so equa- tions 3.39 and 3.40 hold. From 3.43 we find that P1 must satisfy (3.62)∇2P1 = −4pi N ∑ j=1 ( δj1 − P0 ) ajδ ( x − xj ) , x ∈ Ω, ∂nP1 = −λ0P0, x ∈ ∂Ω. The solvability condition for P1 obtained by the divergence theorem gives (3.63)−4pi N ∑ j =1 ( δj1 − P0 ) aj = ∫∫ ∂Ω λ0P0 = −λ0P0|∂Ω|, where |∂Ω| is the surface area of the sphere Ω. Thus P0 = 4pia1 4piNa¯ + λ0|∂Ω| . (3.64) 673.3. Multiple targets with asymptotically null volume Next we decompose P1 = −λ0P0P1P +P1H + χR where χR is a constant and P1P and P1H obey the following conditions ∇2P1P = |∂Ω| |Ω| , x ∈ Ω, (3.65)∂nP1P = 1, x ∈ ∂Ω,∫ Ω P1P (x) dx = 0; ∇2P1H = λ0P0|∂Ω| |Ω| − 4pi N ∑ j=1 ( δj1 − P0 ) ajδ ( x − xj ) , x ∈ Ω, (3.66)∂nP1H = 0, x ∈ ∂Ω,∫ Ω 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 P1H = 4pi N ∑ j=1 ( δj1 −P0 ) ajG ( x; xj ) , (3.67) by identity 3.63. As for P1P, the divergence theorem yields∫ Ω ( G (ξ; x)∇2P1P (ξ) − P1P (ξ)∇2G (ξ; x) ) dξ = ∫ ∂Ω G (ξ; x) dξ. The LHS is equal to λ0P0|∂Ω| |Ω| ∫ Ω G (ξ; x) dξ − ∫ Ω P1P (ξ) ( 1 |Ω| − δ (ξ − x) ) dξ = P1P (x) . Hence P1P (x) = ∫ ∂Ω G (ξ; x) dξ. (3.68) Indeed, by reciprocity of G the third condition for P1P is satisfied too since∫ Ω P1P (x) dx = ∫ ∂Ω ( ∫ Ω G (ξ; x) dx ) dξ = ∫ ∂Ω ( ∫ Ω G (x; ξ) dx ) dξ = 0. 683.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 |∂Ω| |Ω| = 4pia2 4 3 pia 3 = 3 a . Thus the general solution for P1P is P1P (r) = c1 r + c2 + r2 2a where c1 and c2 are constants. Then the boundary condition yields 1 = P ′1P (a) = ( − c1 r2 + r a ) |r=a= − c1 a2 + 1 so that c1 = 0. On the other hand, the integral condition imposed on P1P gives 0 = ∫ Ω P1P (x) dx = 4pi ∫ a 0 ( c2 + r2 2a ) r2dr = 4pi ( c2a 3 3 + a5 10a ) by a change to spherical coordinates. Then c2 = − 3 10 a and so P1P (x) = |x| 2 2a − 3 10 a. (3.69) We have reached the analogous step to equation 3.51, a two term ex- pansion for the splitting probability: (3.70) P (x) = a1Na¯ + λ0a2 + ea1 [ − λ0 Na¯ + λ0a2 ( |x|2 2a − 3 10 a ) + 4piG (x; x1) − 1 Na¯ + λ0a2 N ∑ j=1 ajG ( x; xj )] + eχR +O ( e2 ) . The derivation of χR has been left for Appendix A. See Figures 3.17 and 3.18 for some example plots of P . 693.3. Multiple targets with asymptotically null volume 0 0.2 0.4 0.6 0.8 1 −0.5 0 0.5 1 1.5 2 2.5 3 x 10 −3 r P 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 x 10 −3 r P Figure 3.17: Splitting probability from equations 3.70, A.11, and A.19, with e = 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). 0 1.5708 3.1416 4.7124 6.2832 0.5 1 1.5 2 2.5 3 3.5 x 10 −3 Azimuthal angle P 0 1.5708 3.1416 4.7124 6.2832 2 3 4 5 6 7 8 9 10 x 10 −4 Azimuthal angle P Figure 3.18: Splitting probability from equations 3.70, A.11, and A.19, with e = 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). 703.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 con- dition by the Robin boundary condition we proposed here ∂nTG = eλ0TG. Because this general MFPT is not what we are interested in, I did not in- clude 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 dif- fusion 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 correspondingAPC, we created a mathematical model where the lymph node is layered in all three dimensions. Diffusion on the x direction (respec- tively 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.) 713.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 (3.71)Le [TeP] = P (x) , with P satisfying the corresponding equation from the previous sections. Here e 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 Te in the homogenization limit e → 0, i.e. when the strips are so thin that every region, no matter how small, intersects both types of layers. Notice that Le can be written as Le = −A (x/e) : D2 where D2 is the matrix of second partial derivatives, i.e. with i, j entry equal to ( D2 ) i,j = ∂2 ∂xixj , and the colon operator acts in the following way M : N = tr ( MTN ) = 3 ∑ i,j=1 Mi,jNi,j 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 (3.72)LeTe = 1 [16], where (3.73)Le = − 3 ∑ i,j=1 Ai,j ( x e ) ∂2 ∂xixj , A being the matrix of diffusion coefficients and e the width of the strips. According to [9], in the limit when e → 0 solutions of the general equa- tion Leue = f (x) converge to solutions of the averaged equation − ¯A : D2u = f (x) (3.74) Here the averaged coefficient matrix ¯A is given by ¯A = 〈ρ∞A〉 where the angle brackets denote the average 〈v〉 = ∫ v, and ρ is the invariant distribu- tion of L, that is the solution to L∗ρ∞ = 0, 〈ρ∞〉 = 1. 723.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 〈v〉 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. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 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 (2piz) + 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 corre- sponding to a color strip. In terms of such a function g we could define the diffusion coefficient matrix to be A (x) = D  g(x2)g(x3 − 0.5) 0 00 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 ex- pressed as a product of diagonal matrices each depending on only one xi: A (x) = D n ∏ j=1 Aj ( xj ) . 733.4. Space dependant diffusivity In the particular case we are considering A1 (x1) =  1 0 00 g(x1 − 0.5) 0 0 0 g(x1 − 0.5)  , A2 (x2) =  g(x2) 0 00 1 0 0 0 g(x2)  , A3 (x3) =  g(x3 − 0.5) 0 00 g(x3 − 0.5) 0 0 0 1  . Hence according to [9] the matrix of coefficients for the averaged equation has the simplified form (3.75)¯A = 〈 1 D∏nj=1 Ajjj ( xj )A 〉〈 1 D∏nj=1 Ajjj ( xj ) 〉 −1 = 〈A〉 . If we use g as in Figure 3.20 then h = ∫ 1 0 g (z) dz = ∫ 1 0 g (z− 0.5) dz = 1.12.1 so that the averaged equation ends up being −Dh2∇2u = f (x) (3.76) Notice that the case will be the same whenever the functions that de- scribe the behaviour along each fiber integrate to the same constant, pro- vided 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. 743.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 wherewe consideredmultiple 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 e. (Compare equations 3.51 and 3.70.) Indeed, P0 for the Robin condition has an extra λ0a2 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 e. Moreover, further exten- sions 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 re- 753.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 exam- ple, 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 corre- sponding 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 pre- sented 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 in- fluence of target localization on the killing rates of cytotoxic T cells. Al- though we considered a different process, namely localization of the ac- tivating 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 in- teresting under the T cell activation problem. Following the FRN model, one could also consider the MFPT to a spe- cific 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 direc- tions 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 763.5. Discussion and future work many other biological problems. 77Bibliography [1] Marc Baje´noff, Nicolas Glaichenhaus, and Ronald N. Germain. Fibrob- lastic 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 Be´nichou, BobMeyer, and Raphae¨l Voituriez. First-passage quantities of brownian motion in a bounded domain with multiple targets: a unified approach. Journal of Physics A, Mathe- matical 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. Mathemati- cal 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 retic- ular 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 Mathe- matical Sciences, 7(4), August 2009. [10] Izrail’ M. Gradshteyn and Iosif M. Ryzhik. Table of Integrals, Series, and Products. Academic Press, 1980. 78Bibliography [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. Ro- bust single-particle tracking in live-cell time-lapse sequences. Na- ture 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 Veteri- nary Research, 12(4), July 1964. [15] Theodore Kolokolnikov, Miche`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 exci- tation 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 edi- tion, 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 Univer- sity Press, 2001. [23] Hannes Risken. The Fokker-Planck Equation. Methods of Solution and Ap- plications. Springer-Verlag, second edition, 1989. [24] Ulrich H. von Andrian and Thorsten R. Mempel. Homing and cel- lular 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. 80Appendix 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 equa- tion 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 (A.1)∇2u (x) = −δ (x − ξ) , in R3. Take a small sphere Ωe of radius e about ξ. Then in a neighborhood of ξ we let r = |x− ξ| and then 0 = ∇2u = urr + 2 r ur, for r > 0, so that u (r) = βr . Now, since ∫ Ωe ∇2u (x) dx = − ∫ Ωe δ (x− ξ) dx = −1, then by the divergence theorem −1 = ∫ Ωe ∇u ·~ndS = 4pi ( r2ur (r)) |r=e, (A.2) again by a change to spherical coordinates as in equation 3.45. Since ur = −βr−2 then equation A.2 gives 4pi ( − e2β e2 ) = −1 81A.1. From the Neumann splitting probability, equation 3.51 which yields β = 14pi . Thus u (x) = 14pi|x − ξ| , and so we can write (A.3)G (x; ξ) ≈ 14pi|x − ξ| + R (x; ξ) , as x → ξ. To find R I followed the Appendix A from reference [15]. We first de- compose G as follows G (x; ξ) = u (x) + |x| 2 + |ξ|2 6|Ω| + 1 4pi φ (x; ξ) , (A.4) with φ an unknown function we want to compute. Notice that in that case ∇2G (x; ξ) = ∇2u (x) + 16|Ω| ( ∂2 ∂r2 + 2 r ∂ ∂ ) (|x|2 + |ξ|2) + 14pi∇2φ (x; ξ) = −δ (x − ξ) + 1|Ω| + 1 4pi∇ 2φ (x; ξ) . Because Ω is a sphere, |Ω| = 4pia3/3. Then ∂nG (x; ξ) = 14pi ( ∂n ( 1 |x − ξ| ) + |x| a3 + ∂nφ ) . Substituting this onto the differential equation which defines G, equation 3.49, we find that φ must satisfy (A.5)∇2φ = 0, x ∈ Ω, ∂nφ = − 1 a2 − ∂n ( 1 |x − ξ| ) , x ∈ ∂Ω. 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=0 BnPn (cos θ) ( |x||ξ| a2 ) n , (A.6) 82A.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 en- forced 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φ|∂Ω= ∞ ∑ n=0 nBn an+1 Pn (cos θ) |ξ|n . (A.7) Now, from the generating function for Legendre polynomials, 1√ 1− 2tz + t2 = ∞ ∑ k=0 tkPk (z) , (A.8) we obtain the following equation, by setting z = cos θ and t = |ξ|/|x|: 1 |x − ξ| = 1√|ξ|2 − 2|x||ξ| cos θ + |x|2 = 1 |x|√t2 − 2tz + 1 = ∞ ∑ n=0 |ξ|n |x|n+1 Pn (cos θ) . By differentiating this expression with respect to r = |x| we get (A.9)∂r ( 1 |x − ξ| ) |r=a = − ∞ ∑ n=0 (n + 1) |ξ|n an+2 Pn (cos θ) . Upon substituting A.7 and A.9 into the boundary condition given in A.5, we obtain ∞ ∑ n=0 ( nBn − n + 1 a ) Pn (cos θ) |ξ| n an+1 = − 1 a2 . Since B0 is arbitrary, we can choose it to be equal to 1/a for convenience. Then the other coefficients must satisfy (A.10)Bn = 1 a + 1 na , for n ≥ 1. Back to the series expansion of φ, equation A.6, and replacing on it the previous expressions for the Bn’s, we find that φ (x; ξ) = 1 a ∞ ∑ n=0 Pn (cos θ) ( |x||ξ| a2 )n + 1 a ∞ ∑ n=1 1 n Pn (cos θ) ( |x||ξ| a2 )n . 83A.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 a ∞ ∑ n =0 Pn (cos θ) ( |x||ξ| a2 )n = a |x| ( a4 |x|2 − 2|ξ| a2 |x| cos θ + |ξ| 2 ) − 1 2 = 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 (β) = ∞ ∑ n=1 1 n Pn (cos θ) βn. 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 ′ (β) = ∞∑ n=1 Pn (cos θ) βn−1 = 1 β ∞ ∑ n=1 Pn (cos θ) βn = 1 β ( ∞ ∑ n=0 βnPn (cos θ)− P0 (cos θ) ) = 1 β ( 1√β2 − 2β cos θ + 1 − 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, I (β) = ∫ β 0 1 s ( 1√ s2 − 2s cos θ + 1 − 1 ) ds = log ( 2 1− β cos θ + √1 + β2 − 2β cos θ ) . We are finally in a position to obtain a final expression for φ and so also for G. φ (x; ξ) = a|x|r′ + 1 a log ( 2a2 a2 − |x||ξ| cos θ + |x|r′ ) , 84A.1. From the Neumann splitting probability, equation 3.51 (A.11) G (x; ξ) = 14pi|x − ξ| + |x|2 + |ξ|2 6|Ω| + a 4pi|x|r′ + 1 4pia log ( 2a2 a2 − |x||ξ| cos θ + |x|r′ ) + B, where B is a constant which will be chosen to satisfy ∫ Ω G = 0. Now we want to show that ∫ Ω G (x; ξ) is independent of ξ and so that we can compute B by looking at the simpler equation ∫ Ω G (x; 0) = 0. We have: 0 = ∫ Ω G ( x; ξ′) [ ∇2G (x; ξ) − 1|Ω| + δ (x − ξ) ] dx = ∫ Ω G ( x; ξ′)∇2G (x; ξ) dx − 1|Ω| ∫ Ω G ( x; ξ′) + G (ξ; ξ′) = ∫ Ω ∇ · [ G ( x; ξ′)∇G (x; ξ)]dx − ∫ Ω [ ∇G (x; ξ) · ∇G (x; ξ′)]dx − 1 |Ω| ∫ Ω G ( x; ξ′) + G (ξ; ξ′) . Then the divergence theorem gives∫ Ω ∇ · [ G ( x; ξ′)∇G (x; ξ)]dx = ∫ ∂Ω G ( x; ξ′)∇G (x; ξ) ·~ndx = 0, by the boundary condition on G, i.e. ∂nG = 0 on ∂Ω. Hence 1 |Ω| ∫ Ω G ( x; ξ′) = G (ξ; ξ′) − ∫ Ω ∇G (x; ξ) · ∇G (x; ξ′) dx. Note that the RHS is symmetric in ξ and ξ′. It follows that ∫Ω G (x; ξ) dx =∫ Ω G (x; ξ′) dx. With a change to spherical coordinates, like in equation 3.45, we compute the desired integral: 0 = ∫ Ω G (x; 0) dx = 1 4pi ∫ 2pi 0 ∫ pi 0 ∫ ∞ 0 r sin θdrdθdϕ + 14pia ∫ 2pi 0 ∫ pi 0 ∫ ∞ 0 r2 sin θdrdθdϕ + 1 6|Ω| ∫ 2pi 0 ∫ pi 0 ∫ ∞ 0 r4 sin θdrdθdϕ + B|Ω| = ∫ a 0 rdr + 1 a ∫ a 0 r2dr + 4pi6|Ω| ∫ a 0 r4dr + B|Ω| = a2 ( 1 2 + 1 3 + 1 10 ) + 4pia3 3 B. 85A.1. From the Neumann splitting probability, equation 3.51 Here we have used the fact that when ξ = 0, r′ = |x′| = a2/|x|. Thus B = − 710pia . (A.12) A.1.2 Derivation of χ1 In this section we compute the constant from the order e function, P1, from the splitting probability with multiple targets and Neumann outer bound- ary. Let us first write P1 (x) as x → xj using what we found on the previous section (equation A.3). P1 (x) ≈ (1− P0) a1|x − x1| + 4pi (1− P0) a1R1,1 − 4piP0 N ∑ i=2 aiG1,i + χ1, as x → x1, P1 (x) ≈ − P0aj|x − xj| − 4piP0ajRj,j + 4pia1Gj,1 − 4piP0 N ∑ i=1,i 6=j aiGj,i + χ1, as x → xj, j = 2, . . . , N, where Gi,j ≡ G ( xi; xj ) and Ri,j ≡ R ( xi; xj ) . We reduce the previous equa- tions by introducing the notation B1 = 4pia1R1,1 − 4piP0 ( a1R1,1 + N ∑ i=2 aiG1,i ) , Bj = 4pia1Gj,1 − 4piP0 ( ajRj,j + N ∑ i=1,i 6=j aiGj,i ) , j = 2, . . . , N. Thus P1 (x) ≈    (1−P0) a1 |x− x1| + B1 + χ1, as x → x1 − P0aj |x− x1| + Bj + χ1, as x → xj . (A.13) 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 w1 (y) = c ( 1− aj |y| ) , where c is a constant which is readily found to be c = Bj + χ1, since w1 → c as |y| → ∞. 86A.1. From the Neumann splitting probability, equation 3.51 Going back again to the matching condition 3.41, it yields P2 ≈ − aj ( Bj + χ1 ) |x − xj| , as x → xj, j = 1, . . . , N. By an argument like that led by equations 3.44 and 3.45 we find that P2 should satisfy (A.14)∇2P2 = 4pi N ∑ j=1 aj ( Bj + χ1 ) δ ( x − xj ) , x ∈ Ω, ∂nP2 = 0, x ∈ ∂Ω. The divergence theorem yields a solvability condition for P2, just as it did for P1 on equation 3.47, namely N ∑ j=1 aj ( Bj + χ1 ) = 0. Hence χ1 = − 1 Na¯ N ∑ j=1 ajBj. The expressions for Bj can be substituted back and so χ1 can be written in matrix form as χ1 = − 4pia1 Na¯ [( GTa ) 1 − 1 Na¯ a TGa ] , (A.15) where a = (a1, . . . , aN)T and G ≡   R1,1 G1,2 · · · G1,N G2,1 . . . . . . . . . . . . . . . . . . GN−1,N GN,1 · · · GN,N−1 RN,N   . Notice that with equations A.11 and A.12 we can compute all the entries of this matrix. In particular (A.16)Ri,i = |xi| 2 4pia3 + a 4pi ( a2 − |xi|2) + 1 4pia log ( a2 a2 − |xi|2 ) − 7 10pia . 87A.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 e function, P1, from the splitting probability with multiple targets and Robin outer boundary. First we look at the limit as x → xj of P1 (x). Like in Section A.1.2 we write P1 (x) ≈    (1−P0) a1 |x− x1| + B1 + χR, as x → x1 − P0aj |x− x1| + Bj + χR, as x → xj (A.17) with B1 = −λ0P0P1P (x1) + 4pia1R1,1 − 4piP0 ( a1R1,1 + N ∑ i=2 aiG1,i ) , Bj =−λ0P0P1P ( xj ) +4pia1Gj,1− 4piP0 ( ajRj,j + N ∑ i=1,i 6=j aiGj,i ) , j = 2, . . . , N. Following an argument entirely analogous to that on A.1.2 we find that P2 must satisfy (A.18)∇2P2 = 4pi N ∑ j=1 aj ( Bj + χR ) δ ( x − xj ) , x ∈ Ω, ∂nP2 = −λ0P1, x ∈ ∂Ω. Upon using the divergence theorem on the previous equation we obtain a solvability condition for P2, namely 4pi N ∑ j =1 aj ( Bj + χR ) = ∫ Ω ∇2P2dV = −λ0 ∫ ∂Ω P1dS = λ20P0 ∫ ∂Ω P1PdS − λ04pi N ∑ j=1 ( δj1 − P0 ) ajP1P ( xj ) − λ0χR|∂Ω| 88A.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∫ ∂Ω P1PdS = a|∂Ω| 5 = 4pia3 5 . Substituting back the expressions for Bj we find that∑Nj=1 ajBj can bewritten with matrix notation as N ∑ j =1 ajBj = −λoP0 N ∑ j=1 ajP1P ( xj ) + 4pi [ a1 ( GTa ) 1 − P0aTGa ] . Recall that P1P (x) = |x| 2 2a − 3 10 a. Therefore χR = 1 4pi ( Na¯ + λ0a2 ) ( λ20P0 ∫ ∂Ω P1PdS − λ04pi N ∑ j=1 ( δj1 − P0 ) ajP1P ( xj ) − 4pi N ∑ j=1 ajBj ) = 1 Na¯ + λ0a2 ( 1 5λ 2 0P0a3 − λ0a1P1P (x1) + λ0P0 N ∑ j=1 ajP1P ( xj ) + λoP0 N ∑ j=1 ajP1P ( xj ) − 4pi [ a1 ( GTa ) 1 − P0aTGa ]) = 1 Na¯ + λ0a2 ( λ20P0a3 5 − λ0a1|x1|2 2a + 3λ0a1a 10 + λ0P0 a N ∑ j=1 aj|xj|2 − 3λ0P0aNa¯ 5 − 4pi [ a1 ( GTa ) 1 − P0aTGa ]) . (A.19) 89Appendix B Programs In this appendix I present the most relevant codes written by me which I used for Chapter 2. 90B .1 . sp otD etecto r3D B.1 spotDetector3D Extension to three dimensions of the program spotDetector, algorithm from [21], coded by Franc¸ois Aguet, which can be downloaded from http://lccb.hms.harvard.edu/software.html. 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 91B .1 . sp otD etecto r3D % 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. if nargin<5 dthreshold = 5; end if nargin<4 choice=0; elseif isempty(choice) choice=0; end92B .1 . sp otD etecto r3D 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; 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 %% 3D process connected components mask3d=img3d>0; localMax = locmax3d(img3d, 7); [labels, nComp] = bwlabeln(mask3d, 6); area = zeros(nComp, 1);93B .1 . sp otD etecto r3D 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); xmax2 = cell(nComp, 1); ymax2 = cell(nComp, 1); zmax2 = cell(nComp, 1); area2 = cell(nComp, 1); totalInt2 = cell(nComp, 1); labelVect2 = cell(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); 94B .1 . sp otD etecto r3D 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); values = img3d(stats(n).PixelIdxList); totalInt(n) = sum(values); 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;95B .1 . sp otD etecto r3D 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)); % remove highest max from list xm(maxIdx(1)) = []; ym(maxIdx(1)) = []; zm(maxIdx(1)) = []; % 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); % 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);96B .1 . sp otD etecto r3D 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})); %intensity values totalInt2{n} = totalInt(n)*ones(nSecMax,1)/nMaxima(n); totalInt(n) = totalInt(n)/nMaxima(n); end end end xmax2 = vertcat(xmax2{:}); ymax2 = vertcat(ymax2{:}); zmax2 = vertcat(zmax2{:}); totalInt2 = vertcat(totalInt2{:}); area2 = vertcat(area2{:}); labelVect2 = vertcat(labelVect2{:}); % assign results to output structure frameInfo.xmax = [xmax; xmax2(:)]; frameInfo.ymax = [ymax; ymax2(:)]; frameInfo.zmax = [zmax; zmax2(:)]; frameInfo.xcom = [xcom; xmax2(:)];97B .1 . sp otD etecto r3D 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; % 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); 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;98B .2 . m ySp otD etecto r end frameInfo.path = []; frameInfo.maskPath = []; end parent=fileparts(pwd); if choice==0 sfx=’max’; elseif choice==1 sfx=’com’; end save([parent,filesep,’Results’,filesep,base,’sd3D_’,sfx,’.mat’],’movieInfo’); B.2 mySpotDetector Alternative 3D detection algorithm of my own autorship. 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 regions99B .2 . m ySp otD etecto r % 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.100B .2 . m ySp otD etecto r 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 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 fprintf(’Starting analysis of time point %d...\n’,t); % 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;101B .2 . m ySp otD etecto r 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 if ((bb(4)==1) || (bb(5)==1) ||(bb(6)==1)) % if the region is really of dim less than 3, then there’s no % problem determining its center 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 flag=0; %this flag will turn on if in some plane of the current %3D region there’s more than one 2D connected region 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 % regionprops102B .2 . m ySp otD etecto r % 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 % 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;103B .2 . m ySp otD etecto r 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; 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 current104B .2 . m ySp otD etecto r % plane which had empty intersections with all previous % regions, flaglist helps here %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 if flagj==0 PIL{m}=prevPIL{j}; plane{m}=prevplane{j}; m=m+1; end end % keep 2D regions in current plane that matched no % previous regions105B .2 . m ySp otD etecto r for k=flaglist PIL{m}=newlist{k}; plane{m}=z; m=m+1; end prevPIL=PIL; prevplane=plane; end % end of the current plane of a given 3D region % 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 106B .2 . m ySp otD etecto r 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 PIL={}; plane={}; lidx={}; end %ends case more than one z-slice end % end of current 3D region end %end of current time image parent=fileparts(pwd); if bt==0 sfx=’bup’; 107B .3 . d etecttr 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) if nargin<4 d=0; end prnt=pwd; [~,name]=fileparts(path); if d==0 mkdir(path,’lessnoise’) end movieInfo=repmat(struct(’xCoord’,[],’yCoord’,[],’zCoord’,[],’amp’,[],’height’,[]),T,1); H=[]; for t=1:T cd([prnt,filesep,’SpotDetector’]) tstr=num2str(t,’%3.3d’); 108B .3 . d etecttr 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):...109B .3 . d etecttr 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 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 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’); Function scriptTrackGeneral is called in the following way: 110B .3 . d etecttr 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 column111B .4 . testd etect %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 B.4 testdetect 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 112B .4 . testd etect % 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)); %% 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); %% 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 sphere113B .4 . testd etect 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 end %rescale centers C(:,3)=(C(:,3)-r)/4+1; % add noise to original image % noise=uint8(64*randn(size(imagen))); % imagen=imagen+noise;114B .4 . testd etect %% Spot Detector -- Image denoising % for i=1:nslices % detectionMask = denoise(imagen(:,:,i)); % end %% 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); %% Sort results [~,ord2]=sort(C(:,3));C=C(ord2,:);115B .5 . testtrack [~,ord1]=sort(output(:,3));output=output(ord1,:); %% 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; cd /home/monica/Thesisv2/simulations % this function needs to be improved in the following ways: % - find a better way to quantify error % - check the scaling when measuring error, I think it would make more % sense to rescale movieInfo instead of the centers, but one needs to be % careful with the offsets (ie pixels start at 0.5, not 0) B.5 testtrack 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 116B .5 . testtrack % for no_of_steps seconds. curr=pwd; [prnt]=fileparts(curr); movieInfo=repmat(struct(’xCoord’,[],’yCoord’,[],’zCoord’,[],’amp’,[]),no_of_steps,1); %% Specify simulation parameters dim = 3; % Dimensionality tau = 6; % Sampling interval in s (0.001 = 1 ms) x= 250; z= 24; %% 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); % 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));117B .5 . testtrack % 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)]; % end % 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 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) 118Appendix C Additional plots 0 Pi/4 Pi/2 3Pi/4 Pi0 1000 2000 3000 4000 5000 LN2sc − Turning angle Angle Frequenc y 0 5 10 15 20 25 0 2000 4000 6000 8000 10000 12000 14000 16000 LN2sc − Step size Distance (um) Frequenc y 0 2 4 6 8 10 12 0 50 100 150 200 250 300 350 LN2sc − Average step size Distance (um) Frequenc y 0 50 100 150 0 500 1000 1500 2000 2500 3000 3500 LN2sc − Diameter Distance (um) Frequenc y 0 10 20 30 40 50 60 0 500 1000 1500 2000 LN2sc − Track length Time steps Frequenc y Figure C.1: Histograms for data set LN2 using mySpotDetector-bup as de- tection method with the lower level of denoising. 119Appendix C. Additional plots 0 Pi/4 Pi/2 3Pi/4 Pi0 500 1000 1500 2000 LN8 − Turning angle Angle Frequenc y 0 5 10 15 20 0 2000 4000 6000 8000 10000 12000 14000 LN8 − Step size Distance (pixels) Frequenc y 0 2 4 6 8 10 0 50 100 150 200 250 300 350 LN8 − Average step size Distance (pixels) Frequenc y 0 20 40 60 80 0 500 1000 1500 2000 2500 3000 LN8 − Diameter Distance (pixels) Frequenc y 0 10 20 30 40 50 60 0 500 1000 1500 2000 2500 3000 LN8 − Track length Time steps Frequenc y Figure C.2: Histograms for data set LN8 using mySpotDetector-bup as de- tection method with the lower level of denoising. 120Appendix C. Additional plots Pi/4 Pi/2 3Pi/4 Pi0 200 400 600 800 1000 1200 LN9 − Turning angle Angle Frequenc y 0 5 10 15 20 0 1000 2000 3000 4000 5000 6000 7000 LN9 − Step size Distance (pixels) Frequenc y 0 2 4 6 8 10 12 0 50 100 150 200 250 300 LN9 − Average step size Distance (pixels) Frequenc y 0 10 20 30 40 50 60 0 500 1000 1500 2000 LN9 − Diameter Distance (pixels) Frequenc y 0 5 10 15 20 0 200 400 600 800 1000 1200 1400 LN9 − Track length Time steps Frequenc y Figure C.3: Histograms for data set LN9 using mySpotDetector-bup as de- tection method with the lower level of denoising. 121Appendix C. Additional plots Pi/4 Pi/2 3Pi/4 Pi0 200 400 600 800 z1sc − Turning angle Angle Frequenc y 0 5 10 15 20 0 1000 2000 3000 4000 5000 z1sc − Step size Distance (um) Frequenc y 0 2 4 6 8 10 0 50 100 150 200 250 z1sc − Average step size Distance (um) Frequenc y 0 10 20 30 40 50 0 200 400 600 800 1000 z1sc − Diameter Distance (um) Frequenc y 0 10 20 30 40 0 500 1000 1500 2000 2500 z1sc − Track length Time steps Frequenc y Figure C.4: Histograms for data set z1 using mySpotDetector-bup as detec- tion method with the lower level of denoising. 122Appendix C. Additional plots Pi/4 Pi/2 3Pi/40 5 10 15 20 25 LN10 detection with utrack − angles Angle Frequenc y 0 2 4 6 8 10 0 50 100 150 200 LN10 detection with utrack − step size Distance (um) Frequenc y 0 5 10 15 20 25 0 20 40 60 80 100 LN10 detection with utrack − diameter Distance (um) Frequenc y 0 5 10 15 20 25 30 0 20 40 60 80 100 120 LN10 detection with utrack − length Time steps Frequenc y 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. 123Appendix C. Additional plots Pi/4 Pi/2 3Pi/40 5 10 15 20 25 30 LN10 detection with utrack − angles Angle Frequenc y 0 2 4 6 8 10 0 50 100 150 200 LN10 detection with utrack − step size Distance (um) Frequenc y 0 5 10 15 20 25 0 20 40 60 80 100 120 LN10 detection with utrack − diameter Distance (um) Frequenc y 0 5 10 15 20 25 30 0 20 40 60 80 100 120 140 LN10 detection with utrack − length Time steps Frequenc y 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. 124Appendix C. Additional plots Pi/4 Pi/2 3Pi/40 5 10 15 20 25 30 LN10 detection with utrack − angles Angle Frequenc y 0 2 4 6 8 10 12 0 50 100 150 200 250 LN10 detection with utrack − step size Distance (um) Frequenc y 0 5 10 15 20 25 0 20 40 60 80 100 120 LN10 detection with utrack − diameter Distance (um) Frequenc y 0 5 10 15 20 25 30 0 20 40 60 80 100 120 140 LN10 detection with utrack − length Time steps Frequenc y 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 pos- sibility of immediate direction reversal. 125Appendix C. Additional plots Pi/4 Pi/2 3Pi/40 50 100 150 200 250 300 LN2trd − Turning angle Angle Frequenc y 0 5 10 15 20 0 500 1000 1500 2000 2500 3000 3500 4000 LN2trd − Step size Distance (um) Frequenc y 0 10 20 30 40 50 0 200 400 600 800 1000 1200 LN2trd − Diameter Distance (um) Frequenc y 0 20 40 60 80 100 120 0 500 1000 1500 2000 2500 LN2trd − Track length Time steps Frequenc y 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. 126Appendix C. Additional plots Pi/4 Pi/2 3Pi/40 10 20 30 40 50 LN8trd − Turning angle Angle Frequenc y 0 5 10 15 20 0 100 200 300 400 500 600 700 LN8trd − Step size Distance (um) Frequenc y 0 5 10 15 20 25 0 50 100 150 200 250 LN8trd − Diameter Distance (um) Frequenc y 0 10 20 30 40 50 60 0 50 100 150 200 250 300 350 400 LN8trd − Track length Time steps Frequenc y 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. 127Appendix C. Additional plots Pi/4 Pi/2 3Pi/40 5 10 15 LN9trd − Turning angle Angle Frequenc y 0 2 4 6 8 10 0 50 100 150 200 LN9trd − Step size Distance (um) Frequenc y 0 5 10 15 20 0 10 20 30 40 50 60 70 80 LN9trd − Diameter Distance (um) Frequenc y 0 5 10 15 20 0 20 40 60 80 100 120 140 160 LN9trd − Track length Time steps Frequenc y Figure C.10: These histograms correspond to the analysis of the endoge- nous data set LN9, using u-track as method of detection and with the higher level of denoising. 128Appendix C. Additional plots Pi/4 Pi/2 3Pi/40 1 2 3 4 5 z1trd − Turning angle Angle Frequenc y 0 10 20 30 40 0 10 20 30 40 50 z1trd − Step size Distance (um) Frequenc y 0 10 20 30 40 0 5 10 15 20 z1trd − Diameter Distance (um) Frequenc y 0 5 10 15 20 25 0 10 20 30 40 50 z1trd − Track length Time steps Frequenc y Figure C.11: These histograms correspond to the analysis of the endoge- nous data set Z1, using u-track as method of detection and with the higher level of denoising. 129

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
United States 16 0
China 11 0
Russia 9 0
Canada 3 0
Poland 3 0
Germany 2 13
Ukraine 2 0
Chile 1 0
Romania 1 0
United Kingdom 1 0
Japan 1 0
City Views Downloads
Unknown 22 13
Shenzhen 6 0
Beijing 5 0
Ashburn 4 0
Surrey 3 0
Washington 2 0
Saint Petersburg 2 0
Los Angeles 1 0
Fort Worth 1 0
Santiago 1 0
London 1 0
Sidra 1 0
Tokyo 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0073088/manifest

Comment

Related Items