Mathematical Modelling of Partially AbsorbingBoundaries in Biological SystemsbySarafa Adewale IyaniwuraB.Sc., University of Ilorin, 2012M.Sc., AIMS-Stellenbosch University, 2014A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Mathematics)The University of British Columbia(Vancouver)August 2016c© Sarafa Adewale Iyaniwura, 2016AbstractThis project presents a mathematical framework for identifying partially permeablebiological boundaries, and estimating the rate of absorption of diffusing objects atsuch a boundary based on limited experimental data. We used partial differentialequations (PDEs) to derive probability distribution functions for finding a particleperforming Brownian motion in a region. These distribution functions can be fitto data to infer the existence of a boundary. We also used the probability distri-bution functions together with maximum likelihood estimation to estimate the rateof absorption of objects at the boundaries, based on simulated data. Furthermore,we consider a switching boundary and provide a technique for approximating theboundary with a partially permeable boundary.iiPrefaceThis thesis is original, unpublished, independent work by the author, S. A. Iyani-wura under the supervision of Professor Daniel Coombs.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Biological boundaries . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Definition of terms . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.1 Single Particle Tracking (SPT) . . . . . . . . . . . . . . . 21.2.2 Macrophage . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.3 CD45 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2.4 Immunoglobulin G (IgG) . . . . . . . . . . . . . . . . . . 31.3 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.4 Brief description of chapters . . . . . . . . . . . . . . . . . . . . 62 Mathematical Problem Formulation . . . . . . . . . . . . . . . . . . 82.1 Random walk and the diffusion equation . . . . . . . . . . . . . . 82.2 Mathematical formulation . . . . . . . . . . . . . . . . . . . . . 102.3 Partially absorbing boundary and Robin boundary condition . . . 112.4 Maximum Likelihood Estimation . . . . . . . . . . . . . . . . . . 12iv3 Identifying Biological Boundaries . . . . . . . . . . . . . . . . . . . 133.1 Unbounded domain: Imaginary boundary . . . . . . . . . . . . . 133.2 Bounded domain: Perfectly absorbing boundary . . . . . . . . . . 203.2.1 The half-plane problem . . . . . . . . . . . . . . . . . . . 213.2.2 The disk problem . . . . . . . . . . . . . . . . . . . . . . 264 Partially Absorbing Boundary . . . . . . . . . . . . . . . . . . . . . 324.1 Half-plane problem . . . . . . . . . . . . . . . . . . . . . . . . . 324.2 The disk problem . . . . . . . . . . . . . . . . . . . . . . . . . . 434.2.1 Estimating the rate of absorption on the boundary . . . . . 515 Fluctuating Boundary . . . . . . . . . . . . . . . . . . . . . . . . . . 575.1 Estimating the rate of absorption on a fluctuating boundary . . . . 575.2 Switching boundary . . . . . . . . . . . . . . . . . . . . . . . . . 625.2.1 Switching between perfectly reflecting and partially ab-sorbing . . . . . . . . . . . . . . . . . . . . . . . . . . . 625.2.2 Switching between perfectly reflecting and perfectly ab-sorbing . . . . . . . . . . . . . . . . . . . . . . . . . . . 686 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . 726.1 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78A Supporting Materials . . . . . . . . . . . . . . . . . . . . . . . . . . 81A.1 Smoldyn code for partially absorbing boundary problem . . . . . 81A.2 Smoldyn code for fluctuating boundary problem . . . . . . . . . . 83vList of FiguresFigure 1.1 An example of a trajectory from single particle tracking exper-iment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3Figure 1.2 An illustration of the patterned coverslip . . . . . . . . . . . . 4Figure 1.3 A sketch of one of the circular patterns on the coverslip and adiffusing molecule. . . . . . . . . . . . . . . . . . . . . . . . 5Figure 2.1 A random walk on one dimensional lattice . . . . . . . . . . . 9Figure 3.1 A particle performing Brownian motion in an unbounded two-dimensional domain. . . . . . . . . . . . . . . . . . . . . . . 14Figure 3.2 A particle performing Brownian motion in the Cartesian planewith an imaginary boundary at {(x,y)|x = a,−∞< y< ∞}. . . 16Figure 3.3 The probability of finding a diffusing particle in one of the twohalf-planes of the Cartesian plane at time t . . . . . . . . . . 17Figure 3.4 A particle performing Brownian motion in an unbounded do-main, with an imaginary boundary bounding a disk-shaped re-gion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18Figure 3.5 A sketch of the disk-shaped region. . . . . . . . . . . . . . . 18Figure 3.6 The probability of finding a particle performing Brownian mo-tion inside a disk-shaped region of radius a with an imaginaryboundary at time t . . . . . . . . . . . . . . . . . . . . . . . 20Figure 3.7 A particle performing Brownian motion in the right half-planewith an absorbing boundary at x = 0. . . . . . . . . . . . . . 21Figure 3.8 An illustration of the idea of method of images. . . . . . . . 22viFigure 3.9 The probability of finding a particle performing a Brownianmotion in the right-half plane at time t . . . . . . . . . . . . . 25Figure 3.10 A particle performing Brownian motion in a disk-shaped re-gion with a perfectly absorbing boundary. . . . . . . . . . . . 26Figure 3.11 A sketch of the method of images for a disk. . . . . . . . . . 28Figure 3.12 The probability of finding a particle performing Brownian mo-tion inside a disk-shaped of radius a at time t . . . . . . . . . 30Figure 4.1 A particle performing Brownian motion in the right half-planewith a partially absorbing boundary at x = 0. . . . . . . . . . 33Figure 4.2 The comparison of the probability of finding the particle in theright half-plane obtained from the simulation and that of theprobability distribution function in Equation (4.11). . . . . . . 40Figure 4.3 The mean and variance of first passage time for a particle per-forming Brownian motion in the right half-plane, with a par-tially absorbing boundary at x = 0. . . . . . . . . . . . . . . . 42Figure 4.4 The comparison of the probability of finding a particle in thedisk-shaped region calculated from simulation and that of theprobability distribution function in Equation (4.36). . . . . . . 49Figure 4.5 The mean and variance of first passage time for a particle per-forming Brownian motion in a disk-shaped region with a par-tially absorbing boundary. . . . . . . . . . . . . . . . . . . . 50Figure 4.6 Estimates of the rate of absorption on the boundary with stan-dard deviation error bars based on simulated data consideringdifferent particle initial positions x0, for the half-plane problem. 53Figure 4.7 Estimates of the rate of absorption on the boundary with stan-dard deviation error bars based on simulated data using differ-ent κ values, for the half-plane problem. . . . . . . . . . . . . 54Figure 4.8 Estimates of the rate of absorption on the boundary with stan-dard deviation error bars based on simulated data consideringdifferent particle initial positions r0, for the disk problem. . . 55viiFigure 4.9 Estimates of the rate of absorption on the boundary with stan-dard deviation error bars based on simulated data using differ-ent κ values, for the disk problem. . . . . . . . . . . . . . . . 55Figure 5.1 A particle performing Brownian motion in some region boundedby some other diffusing particles. . . . . . . . . . . . . . . . 58Figure 5.2 Comparison of the probability of finding the particle of interestin the right half-plane calculated from the simulation and theprediction from the mean field approximation. . . . . . . . . . 61Figure 5.3 Comparison of the probability of finding the particle of interestin the a disk-shaped region calculated from the simulation andthe prediction from the mean field approximation. . . . . . . . 61Figure 5.4 Comparison of the probability of finding the particle of interestin the right half-plane calculated from simulation and that ofthe probability distribution function in Equation (5.7). . . . . 64Figure 5.5 Fitting the probability distribution function in Equation (5.7) toresult calculated from simulation in order to get the effectiverate of absorption on the boundary. . . . . . . . . . . . . . . . 65Figure 5.6 Estimates of the effective rate of absorption κe f f obtained byfitting probability distribution functions to the simulated re-sults for different values of κ2 and switching time. . . . . . . 66Figure 5.7 Estimates of the values of the functions f (α) and g(α) ob-tained by fitting the model in Equation (5.9) to the curves inFigure 5.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . 66Figure 5.8 Fittings functions to the estimated values for the functions f (α)and g(α). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67Figure 5.9 Fitting the probability distribution function in Equation (5.7)to the probability calculated from simulation. . . . . . . . . . 68Figure 5.10 Fitting the probability distribution function in Equation (5.8)to the probability calculated from simulation. . . . . . . . . . 69Figure 5.11 Estimates of the effective rate of absorption κe f f obtained byfitting probability distribution functions to the probability cal-culated from simulation. . . . . . . . . . . . . . . . . . . . . 69viiiFigure 5.12 Fitting the model in Equation (5.13) to the estimated effectiverates of absorption in Figure 5.12. . . . . . . . . . . . . . . . 70ixAcknowledgmentsI would like to express my sincere gratitude to my supervisor Daniel Coombs forbelieving in me, and also for his valuable and constructive suggestions throughoutthe planning and development of this thesis. I wish to acknowledge with muchappreciation his guidance and encouragement during my studies.A special gratitude to Michael Ward and Leah Keshet for taking their time togo through this thesis, and also for their constructive comments and suggestions. Iappreciate your kindness.In addition, I would like to say thank you to all the students, staff, and facultymembers of the Institute of Applied Mathematics and the department of mathe-matics for their support, and also for being a great source of motivation. Thankyou Dan for giving me the opportunity to be part of this wonderful department,institute, and the University of British Columbia in general.Finally, I would like to thank my family and friends for their unlimited supportand encouragement. I am very grateful and appreciate your love. I wish to dedicatethis project to my late parents and my siblings.xChapter 1Introduction1.1 Biological boundariesBiological boundaries are barriers that separate a medium into different regions inthe sense that a substance crossing from one region of the medium to another mayexperience some differences in its environment, and as a result react differently.In cell biology, it is known that proteins such as receptors are constantly diffusingon the cell surface. However, it is sometimes observed in experiments that somereceptors experience restricted motion which leads to the speculation that there isa barrier preventing the free diffusion of these molecules. These barriers can eitherbe formed by some other molecules diffusing on the cell surface, or by moleculesthat are confined to some region of the cell surface. Identifying these boundariesin experimental data is important because it will give us a profound understandingof the motion of the molecules involved. However, this may be challenging dueto the size and quantity of the molecules, and the type of motion they undergo.Despite the importance of identifying these boundaries, not much has been done indeveloping mathematical techniques for accomplishing this task.In this project, we present a mathematical framework for identifying partiallypermeable biological boundaries, and estimating the rate of absorption of diffusingobjects at such a boundary based on limited experimental data. This technique isbased on deriving a probability distribution function for finding a diffusing particlein a region at some specified time, and fitting the distribution function to data in1order to infer the existence of a boundary. The derived distribution function is alsoused together with maximum likelihood estimation to estimate the rate of absorp-tion on the boundary. In addition, we looked at switching boundaries and present atechnique for approximating the boundaries with partially permeable boundaries.1.2 Definition of termsBelow are brief definitions of some of the biological terms used in this thesis.1.2.1 Single Particle Tracking (SPT)Single Particle Tracking (SPT) is a technique used to study the motion of an ob-ject in a medium [12]. This technique involves labelling the particle of interestwith fluorescent optical labels. The sample containing this particle is placed underthe microscope where it is excited using a laser light. When the fluorophores areexcited, several pictures of the sample are taken over a period of time, and eachpicture frame is analyzed to obtain the (x,y) coordinates of the particle of interestat a particular time. These coordinates are then put together in order to create thetrajectory of motion of the particle [9]. The tracks from SPT can be analyzed todetermine the type of motion, and also heterogeneity in the motion of the particle,which gives information about the particle’s interaction with its surrounding. Fig-ure 1.1 shows an example of a trajectory from single particle tracking experiment.1.2.2 MacrophageA macrophage is a type of white blood cell that uses a process called phagocy-tosis to engulf and digest cellular debris, cancer cells, pathogens, and any otherforeign substance that does not have the kind of protein specific of healthy cellson its surface [16]. Macrophages also help in regulating immune responses andparticipate in the development of inflammation by producing chemical substancessuch as enzymes, complement proteins, and interleukin-1, among others [16].2Figure 1.1: An example of a trajectory from single particle tracking experi-ment.1.2.3 CD45CD45 is a transmembrane protein tyrosine phosphatase that is abundantly expressedon the surface of various types of blood cells, including macrophages [11]. It playsan important role in the functioning of these cells, and can act as both positive andnegative regulator of cellular signals [13].1.2.4 Immunoglobulin G (IgG)Immunoglobulin G is the most abundant type of antibody found in all body fluids[14]. It is produced and released by plasma B cells, and protects the body againstbacterial, viral, and fungi infections by binding to the pathogens. An IgG moleculehas a protein complex that is made up of four peptide-chains, comprising of twoidentical heavy chains and two identical light chains arranged in a Y-shape, andhas two binding sites [14].1.3 MotivationThis project is motivated by a Single Particle Tracking (SPT) experiment involvingthe tracking of the membrane phosphatase CD45 on macrophages. The experi-3mental set-up involves seeding macrophages that have CD45 molecules on theircell surface labelled with quantum dots on a coverslip. The coverslip is patternedwith Immunoglobulin G (IgG), and the pattern consists of 2µm diameter circleswith centres separated by 6µm, arranged in a regular 2D array as shown in Figure1.2.Figure 1.2: An illustration of the patterned coverslip. The circular regionsare 2µm in diameter and contain the IgG (in blue). These regions areseparated by a distance of 6µm from the centre of each region.The macrophages are allowed to settle on the patterned coverslip, and whencontact is formed between the macrophages and the IgG-patterned surface, thequantum dot-labelled CD45 molecules on the cell membrane are tracked for 10seconds at a frequency of 30 Hz (hertz). During this experiment, it was observedthat the contact of the macrophages with the immobile IgG creates a “zone of CD45exclusion”. While some of the phosphatase molecules that fall outside this zoneappear to diffuse freely but bounce off the edge of the zone, those that fall insidethe zone are observed to have restricted motion.4(a) A diffusing molecule bouncing off theexclusion zone.(b) A diffusing molecule absorbed into theexclusion zone.Figure 1.3: A sketch of one of the circular patterns on the coverslip and adiffusing molecule.Figure 1.3a shows a trajectory bouncing off the boundary of an exclusion zoneon the coverslip, while Figure 1.3b shows a trajectory that is absorbed into thiszone. The apparent bouncing of phosphatase molecules from the zone of exclusionwas observed from three different trajectories. However, there is not enough evi-dence to conclude whether the observations represent an anisotropic barrier exclu-sion effect on the trajectories, or the molecules are undergoing isotropic Brownianmotion. In other words, the question is: Are the trajectories actually bouncing offthe “zone of CD45 exclusion”or we are observing simple Brownian motion?We approach this problem from the point of view of a particle performingBrownian motion in a domain, and derive the probability distribution functionfor finding the particle in some region after a specified time, considering differ-ent scenarios and geometries. Several previous works have studied the motion ofa Brownian particle in a domain, and have derived probability distribution func-tions for the particle’s position after some specified time [3], [5], [6], [7], and [10].Chandrasekhar [3] derived the probability density function for finding a particleperforming random walk at a point after taking some specified number of steps,from the point of view of a Bernoulli distribution. He considered the case of anunbounded domain, and a bounded domain with a reflecting boundary and an ab-sorbing boundary. In this project, we derive probability distribution functions for a5particle performing Brownian motion in an unbounded domain, in a bounded do-main with an absorbing boundary and a partially absorbing boundary using partialdifferential equations.1.4 Brief description of chaptersChapter twoIn chapter two, we derive the diffusion equation as a continuous limit of a randomwalk, and then used the derived equation to set-up our mathematical problem. Wealso present a brief description of a partially absorbing boundary, Robin boundarycondition, and maximum likelihood estimation.Chapter threeIn this chapter, we derive probability distribution functions for finding a particleperforming Brownian motion in a region of an unbounded and a bounded domain.For the unbounded domain, we consider a semi-infinite rectangular region and adisk-shaped region, both with ‘imaginary’ boundaries, while for the bounded do-main, we consider the right-half plane with an absorbing boundary at x = 0, and adisk-shaped region, also with an absorbing boundary. These probability distribu-tion functions can be fit to data to infer the existence of a boundary.Chapter fourIn chapter four, we derive probability distribution functions for finding a particleperforming Brownian motion in a domain that is bounded by a partially absorb-ing boundary. For this case, we also consider the right-half plane with a partiallyabsorbing boundary at x = 0 and a disk-shaped region. We use these probabilitydistribution functions together with maximum likelihood estimation to estimate therate of absorption of particles on the partially absorbing boundaries. In addition,we calculate the mean and variance of first passage time for the particle performingBrownian motion to exit the regions.6Chapter fiveIn this chapter, we consider scenarios where the boundary of the domain in whichthe particle is performing Brownian motion switches from one boundary type toanother over time. First, we consider the case where the boundary fluctuates be-tween a perfectly reflecting boundary and different partially absorbing boundaries,and then use mean field approximation to replace the fluctuating boundary with apartially absorbing boundary. We also consider a case where the boundary switchesbetween two specific boundary types, and provide a technique for approximatingthe switching boundary with a partially absorbing boundary.Our goal in this project is to develop a mathematical technique for identifyingpartially permeable biological boundaries, and to estimate the rate of absorptionon such boundaries. We also want to present a technique for approximating aswitching boundary with a partially permeable boundary.7Chapter 2Mathematical ProblemFormulationIn this chapter, we derive the diffusion equation as a continuous limit of a randomwalk, and use the equation to formulate our mathematical problem. In addition,we present a brief introduction to a partially absorbing boundary, Robin boundarycondition, and maximum likelihood estimation.2.1 Random walk and the diffusion equationA random walk is a path consisting of a succession of independent random steps.For example, the path traced by a molecule as it travels in a liquid or gas, the searchpath of a foraging animal, and the price of a fluctuating stock can all be modeledas random walks, although in reality, they may not be random (cf. [18]). Here, weshall show the relationship between a random walk and the diffusion equation.Consider a particle performing random walk on a one dimensional lattice. Sup-pose the rule of its movement is that at each time step of size ∆t, the particle caneither jump to the left or right a distance of ∆x with equal probability, 1/2.Let P(x, t) be the probability that the particle is at position x at time t. Thenaccording to the rule of movement of the particle, there are only two possibilitiesthat the particle reaches position x at time t+∆t: The particle was at x−∆x at timet and jumped to the right; or the particle was at x+∆x at time t and jumped to the8left as shown in Figure 2.1.Figure 2.1: A random walk on one dimensional lattice. This figure shows thedifferent possibilities of having a particle at position x at time t +∆t.The particle can either be at position x−∆x at time t and take a stepforward or at position x+∆x at time t and take a step backward.Since the next movement of the particle is independent of its present location,the probability that the particle is at position x at time t +∆t given that it was atposition x−∆x at time t is 12 P(x−∆x, t), while the probability that the particle is atposition x at time t+∆t given that it was at x+∆x at time t is 12 P(x+∆x, t). Thus,the probability that the particle is at position x at time t+∆t isP(x, t+∆t) =12P(x−∆x, t)+ 12P(x+∆x, t) . (2.1)Let us Taylor expand each term in this equation:P(x, t+∆t) = P(x, t)+∂∂ tP(x, t)∆t+O(∆t)2 , (2.2)P(x±∆x, t) = P(x, t)± ∂∂xP(x, t)∆x+12∂ 2∂x2P(x, t)(∆x)2+O(∆x)3 . (2.3)Substituting the Taylor expansions into Equation (2.1) and simplifying, we have∂∂ tP(x, t)∆t+O(∆t)2 =12∂ 2∂x2P(x, t)(∆x)2+O(∆x)3 .Dividing through by ∆t gives∂∂ tP(x, t)+O(∆t) =(∆x)22∆t∂ 2∂x2P(x, t)+O((∆x)3∆t).9Next, we take the limit ∆x −→ 0, ∆t −→ 0 in such a way that lim (∆x)22∆t = D exits,that is, D is finite: This special limit is called the diffusion limit. Thus, we have thediffusion equation∂∂ tP(x, t) = D∂ 2∂x2P(x, t) . (2.4)In n dimensions, similar calculations yield∂∂ tP(−→x , t) = ∇2P(−→x , t) , (2.5)where −→x ≡ (x1, . . . ,xn) is an n dimensional vector, D is the diffusion coefficient,and ∇2 ≡(∂ 2∂x21+ · · ·+ ∂ 2∂x2n)is the Laplace operator.This shows that the diffusion equation is a continuous limit of a random walkprocess. Throughout this project, we shall use the diffusion equation in Equation(2.5) to model Brownian motion.2.2 Mathematical formulationConsider a particle performing Brownian motion in a two dimensional domain Ω.Let P(−→x , t|−→x 0, t0) be the probability that the particle is at position −→x at time tgiven that it started at point −→x 0 at time t0. Then P satisfies∂P∂ t= D∇2P , −→x ∈Ω, t > t0 ,P(−→x , t0|−→x 0, t0) = δ (−→x −−→x 0) , −→x ≡ (x,y) .(2.6)If we let Γ(t|−→x 0, t0) be the probability that the particle is in a region R at timet given that it started at point −→x 0 ∈Ω at time t0, thenΓ(t|−→x 0, t0) =∫RP(−→x , t|−→x 0, t0) dA . (2.7)We shall solve the problem in Equation (2.6) in different geometries and withdifferent types of boundary conditions.102.3 Partially absorbing boundary and Robin boundaryconditionA boundary is said to be partially absorbing if when a particle hits it, the particleis either reflected back into the domain or is absorbed by the boundary. Erbanand Chapman [4] showed that a partially absorbing boundary can be accuratelymodelled using the Robin boundary condition;D∂∂nP(−→x , t) = κP(−→x , t) , (2.8)where κ is the rate of absorption on the boundary, D is the diffusion coefficient ofthe particle, and ∂∂n is the outward normal derivative on the boundary.Note that when κ = 0 in Equation (2.8), the boundary condition correspondsto the perfectly reflecting boundary condition, ∂∂n P(−→x , t) = 0, while as κ tends toinfinity, it becomes the perfectly absorbing boundary condition, that is, P(−→x , t) =0. This shows that the Robin boundary condition is a weighted average of thereflecting and absorbing boundary conditions.In the context of particle-based simulation, Erban and Chapman [4] went fur-ther to develop a relation between the probability of absorption upon hitting theboundary for a particle performing Brownian motion in a bounded domain and therate of absorption in the Robin boundary condition (2.8). Their relation is given asκ = 2P√Dpi ∆t, (2.9)where P is the probability of absorption on the boundary, κ is the rate of absorptionin the Robin boundary condition, and ∆t is the time step used in the simulation.Steven Andrews [1] also presented empirical relations between the Robin ab-sorption rate and the probability of absorption on the boundary, and the relationsareκ ′ =1√2piP+0.24761P2+0.00616P3+0.20384P4 , (2.10)P= κ ′√2pi−3.33321κ ′2+3.35669κ ′3−1.52092κ ′4 , (2.11)11where κ ′ = κ√∆t2D .Unlike the relation presented by Erban and Chapman (in Equation (2.9)) whichwas derived analytically, Andrews relations were derived by fitting polynomials tosimulated data [1].2.4 Maximum Likelihood EstimationA likelihood function is defined as a function of the parameters of a model givensome observations or data [15]. Let P(X;θ) be the probability of having the set ofobservations X given a set of parameters θ , then the likelihood function L (θ |X),of having the parameter values θ given the outcomes X is given as [15]L (θ |X)∼ P(X;θ) .Maximum Likelihood Estimation (MLE) is a method of estimating the parametersof a model based on some observations or data. In order words, for a fixed dataand a corresponding model, MLE selects the set of values of the model parametersthat maximizes the likelihood function.12Chapter 3Identifying Biological BoundariesIn this chapter, we consider a particle performing Brownian motion in an un-bounded domain and derive the probability distribution function for finding theparticle at a position after a specified time. Then, we assume that there is a region ofthe unbounded domain that is bounded by an ‘imaginary’ boundary, and derive theprobability distribution function for finding the particle in this region. The bound-ary is called imaginary because the diffusing particle is assumed to cross it fromboth directions without any restriction. In other words, the boundary does not haveany effect on the motion of the particle. For the region bounded by the imaginaryboundary, we consider a semi-infinite rectangular region and a disk-shaped region.In addition, we derive probability distribution functions for a particle performingBrownian motion in domains that are bounded by perfectly absorbing boundaries.3.1 Unbounded domain: Imaginary boundaryConsider a particle performing Brownian motion in an unbounded two-dimensionaldomain as show in Figure 3.4. If we let Ω be the Cartesian plane, the probabilitythat the particle is at position−→x at time t given that it started at position−→x 0 at timet = 0, P(−→x , t|−→x 0,0) satisfies∂P∂ t= D∇2P , −∞< x,y < ∞, t > 0 ,P(−→x ,0) = δ (−→x −−→x 0) .(3.1)13Figure 3.1: A particle performing Brownian motion in an unbounded two-dimensional domain.We define the Fourier Transform (F.T.) of P(x,y, t) asℑ(k1,k2, t) =1(2pi)2∫ ∞−∞∫ ∞−∞P(x,y, t)e−i(k1x+k2y) dxdy , (3.2)and the inverse F.T. asP(x,y, t) =∫ ∞−∞∫ ∞−∞ℑ(k1,k2, t)ei(k1x+k2y) dk1 dk2 . (3.3)Taking the F.T. of the PDE in Equation (3.1), we haveddtℑ(k1,k2, t) =−D(k21 + k22)ℑ(k1,k2, t) .Solving this Ordinary Differential Equation (ODE) givesℑ(k1,k2, t) = Ae−D(k21+k22)t , (3.4)14where A is a constant.Taking the F.T. of the initial condition in Equation (3.1), we haveℑ(k1,k2,0) =1(2pi)2e−i(k1x0+k2y0) . (3.5)Applying this initial condition to the solution in Equation (3.4) gives,ℑ(k1,k2, t) =1(2pi)2e−D(k21+k22)t−i(k1x0+k2y0) . (3.6)Let us take the inverse F.T. of Equation (3.6) using Equation (3.3).P(x,y, t) =∫ ∞−∞∫ ∞−∞1(2pi)2e−D(k21+k22)t−i(k1x0+k2y0) ei(k1x+k2y) dk1 dk2 .Evaluating the integrals and simplifying givesP(x,y, t|x0,y0,0) = 14piDt exp(− 14Dt[(x− x0)2 + (y− y0)2]). (3.7)This is the probability distribution function for finding a particle performing Brow-nian motion in an unbounded two-dimensional domain at position (x,y) at time tgiven that it started at position (x0,y0) at time t = 0. It is important to point outthat this equation is the free space Green’s function for the diffusion equation intwo-dimensions.Suppose that there is an imaginary boundary at x = a such that the Cartesianplane is divided into two regions as shown in Figure 3.2. We are interested in de-riving the probability distribution functions for finding the particle in these regionsat time t. Specifically, if we let Ω1 = {(x,y)|a≤ x< ∞ ,−∞< y< ∞}, we wantto find the probability distribution function for finding the particle in Ω1 at time tgiven that the particle started the Brownian motion at (x0,y0) ∈ Ω1 at time t = 0.To get this distribution function, we need to integrate Equation (3.7) over Ω1. Thatis,Γ(t|x0) =∫Ω1P(x,y, t|x0,y0,0) dA =∫ ∞−∞∫ ∞aP(x,y, t|x0,y0,0) dxdy .15Figure 3.2: A particle performing Brownian motion in the Cartesian planewith an imaginary boundary at {(x,y)|x = a,−∞< y< ∞}.Substituting P(x,y, t) as given in Equation (3.7) into this equation, we haveΓ(t|x0) = 14piDt∫ ∞−∞∫ ∞aexp(− 14Dt[(x− x0)2 + (y− y0)2])dxdy .Evaluating the integrals gives the probability distribution for finding the particle inΩ1 over time,Γ(t|x0) = 12 erfc(a− x0√4Dt), (3.8)where erfc(x)=1 - erf(x) is the complementary error function.Figure 3.3 shows the plots of the probability distribution function in Equation(3.8). We observe from this figure that for each value of a, the probability offinding the particle in Ω1 starts from one and decreases as time goes on. This isbecause the particle started the Brownian motion in Ω1 and as time goes on, thechances of it crossing the imaginary boundary increases which leads to an increasein the probability of finding the particle in Ω \Ω1 as seen in Figure 3.3b. Wealso notice from the plots in this figure that the probability of finding the particle ineach region after a specific time depends on the position of the imaginary boundary.16(a) Probability of finding the particlein Ω1.(b) Probability of finding the particlein Ω\Ω1.Figure 3.3: The probability of finding a diffusing particle in one of the twohalf-planes of the Cartesian plane at time t, given that it started a Brow-nian motion at the point x0 = 5, with diffusion coefficient D = 0.1, andan imaginary boundary at x = a.Since the particle can cross the boundary from both directions without restriction,the probability of finding the particle in either of the two regions separated by theboundary converges to 0.5 after a long time. This can be seen by taking the limitof the probability distribution function in Equation (3.8) as t tends to infinity asshown below,limt−→∞Γ(t|x0) = limt−→∞12erfc(a− x0√4Dt)=12.Now, suppose the imaginary boundary bounds a disk-shaped region as shownin Figure 3.4. Then we need to integrate Equation (3.7) over the disk-shaped regionto get the probability distribution function for finding the particle in this region attime t.Before integrating, let us write the function in polar coordinates. From Equa-tion (3.7), we haveP(−→x , t|−→x 0,0) = 14piDt exp(−|−→x −−→x 0|24Dt). (3.9)Let |−→x 0| = r0, and |−→x | = r as shown in Figure 3.5. Using the cosine rule, the17Figure 3.4: A particle performing Brownian motion in an unbounded domain,with an imaginary boundary bounding a disk-shaped region.Figure 3.5: A sketch of the disk-shaped region.probability distribution function in Equation (3.9) becomesP(r,θ , t) =14piDtexp(− 14Dt(r2+ r20−2rr0 cos(θ0−θ))). (3.10)Next, we integrate P(r,θ , t) over the disk-shaped region to get the probability dis-18tribution function for finding the particle in the region at time t. That is,Γ(t|r0) = 14piDt∫ 2pi0∫ a0exp(− 14Dt(r2+ r20−2rr0 cos(θ0−θ)))r dr dθ ,Γ(t|r0) = 14piDt e(− r204Dt) ∫ a0r e(− r24Dt)[∫ 2pi0exp(14Dt(2rr0 cos(θ0−θ)))dθ]dr.But ∫ 2pi0exp(12Dt(rr0 cos(θ0−θ)))dθ = 2pi I0( r0r2Dt),where I0(z) is the modified Bessel function of the first kind of order zero.Therefore,Γ(t|r0) = 14piDt exp(− r204Dt)[2pi∫ a0r exp(− r24Dt)I0( r0r2Dt)dr].Evaluating the integral in this equation, we haveΓ(t|r0) = 14piDt exp(− r204Dt)[4piDt exp(r204Dt) (1−Q1(r0√2Dt,a√2Dt))].Simplifying, we have the probability distribution function for finding the particlein the disk-shaped region at time t given that it start at a distance r0 away from thecenter of the region,Γ(t|r0) = 1−Q1(r0√2Dt,a√2Dt), (3.11)where Q1 is the Marcum Q-function of order one.The Marcum Q-function of order M QM is defined asQM(a,b) =∫ ∞bx( xa)M−1exp(−x2+a22)IM−1 (ax) dx , (3.12)where IM−1 is the modified Bessel function of the first kind of order M−1 [17].Figure 3.6 shows the plots of the probability distribution function in Equation(3.11). In this figure, we fixed the initial position of the particle r0 = 5, and varied19(a) Probability of finding the particleinside the disk-shaped region.(b) Probability of finding the particle outsidethe disk-shaped region.Figure 3.6: The probability of finding a particle performing Brownian motioninside a disk-shaped region of radius a with an imaginary boundary attime t, given that it started at a distance r0 = 5 away from the center ofthe region and with diffusion coefficient D = 0.1.the radius of the region. We notice from these plots that as the radius of the disk-shaped region increases, the probability that the particle is still in the region aftera specific time increases. This shows that the probability of finding the particle inthe region at some time t > 0 depends on the radius of the region. This is becausethe time required for the particle to reach the boundary increases as the regiongets bigger. Similar to the plots of the distribution function in Equation (3.8), weobserve that the probability of finding the particle inside the disk-shaped regionstarts from one and decreases with respect to time. The reason for this is that theparticle started its motion inside the region and as time goes on, the chances of itcrossing the imaginary boundary increases. In addition, the probability of findingthe particle inside the disk-shaped region tends to zero after a long time. This canearly be seen by taking the limit as t tends to infinity of the probability distributionfunction in Equation (3.11).3.2 Bounded domain: Perfectly absorbing boundaryIn this section, we derive probability distribution functions for a particle perform-ing Brownian motion in a bounded domain with a perfectly absorbing boundary.20Since the boundary is perfectly absorbing, the diffusing particle vanishes immedi-ately when it hits the boundary. We again consider two different geometries, theright half-plane (rectangular coordinates) and a disk-shaped region (polar coordi-nates).3.2.1 The half-plane problemSuppose we have a particle performing Brownian motion in the right half-planewith a perfectly absorbing boundary at x = 0 as shown in Figure 3.7. Let us defineΩ= {(x,y) |0≤ x< ∞, −∞< y< ∞}. We impose a Dirichlet (absorbing) bound-ary condition on the boundary at x = 0 for the problem in Equation (2.6) to get theproblem for the probability of finding the particle at a point −→x at time t given thatit started at −→x 0 at time t = 0, P(−→x , t|−→x 0);Figure 3.7: A particle performing Brownian motion in the right half-planewith an absorbing boundary at x = 0.21∂P∂ t= D∇2P , −→x ∈Ω , t > 0 ,P(−→x ,0) = δ (−→x −−→x 0) ,P = 0 , x = 0 ,−∞< y< ∞ ,P−→ 0 , as x−→ ∞ .(3.13)We shall rewrite this problem using the method of images which involves extendingthe domain Ω by adding its mirror image with respect to a symmetry hyperplane.Then we place a mirror image of the initial condition in the extended domain insuch a way that the boundary condition P = 0 is satisfied at x = 0. The originalFigure 3.8: An illustration of the idea of method of images.initial condition δ (x− x0,y− y0) is at (x0,y0), while the image initial conditionδ (x+ x0,y− y0) is placed at (−x0,y0) as shown in Figure 3.8. So doing, the new22problem is∂P∂ t= D∇2P, −∞< x,y< ∞, t > 0,P =0, on x = 0, −∞< y< ∞,P =δ (x− x0,y− y0)−δ (x+ x0,y− y0), at t = 0 .(3.14)We observe that our new initial condition is the difference of the actual and imageinitial conditions. This is required in order to satisfy the boundary condition P = 0at x = 0.Define the Fourier Transform (F.T.) of P(x,y, t) as,P̂(k1,k2, t) =1(2pi)2∫ ∞−∞∫ ∞−∞P(x,y, t)e−i(k1x+k2y) dxdy , (3.15)and the inverse F.T. asP(x,y, t) =∫ ∞−∞∫ ∞−∞P̂(k1,k2, t)ei(k1x+k2y) dk1 dk2 . (3.16)Taking the F.T. of the PDE in Equation (3.14),1(2pi)2∫ ∞−∞∫ ∞−∞Pt(x,y, t)e−i(k1x+k2y) dxdy =D1(2pi)2∫ ∞−∞∫ ∞−∞∇2P(x,y, t)e−i(k1x+k2y) dxdy ,P̂(k1,k2, t) =−D(k21 + k22) P̂(k1,k2, t) .Solving this equation, we haveP̂(k1,k2, t) = Ae−D(k21+k22)t , (3.17)where A is a constant.Taking the F.T. of the initial condition in Equation (3.14), we haveP̂(k1,k2,0) =1(2pi)2(e−i(k1x0+k2y0)− e−i(k1x0−k2y0)).23Applying this initial condition to the solution in Equation (3.17), we obtainP̂(k1,k2,σ) =1(2pi)2(e−i(k1x0+k2y0)− e−i(k1x0−k2y0))e−D(k21+k22)t .Now, let us take the inverse F.T. of this equation using Equation (3.16),P(x,y, t) =1(2pi)2∫ ∞−∞∫ ∞−∞e−D(k21+k22)t+ i[k1(x−x0)+k2(y−y0)] dk1 dk2+1(2pi)2∫ ∞−∞∫ ∞−∞e−D(k21+k22)t+ i[k1(x+x0)+k2(y−y0)] dk1 dk2 .Evaluating these integrals using the fact that∫ ∞−∞e−Dk21t+ik(z−x) dk =( piDt)1/2exp((z− x)24Dt),we obtain the probability distribution function for finding a particle performingBrownian motion in the right half-plane at position −→x at time t given that it startedfrom the point −→x 0 at t = 0, and that there is an absorbing boundary at x = 0,P(−→x , t|−→x 0) = 14Dpit[exp(− [(x− x0)2 + (y− y0)2]4Dt)−exp(− [(x+ x0)2 + (y− y0)2]4Dt) ], x> 0 .(3.18)To get the probability distribution function for finding the particle in the right half-plane at time t, we need to integrate P(−→x , t|−→x 0) over the region. That is,Γ(t|x0) =∫ ∞−∞∫ ∞0P(x,y, t|x0,y0) dxdy .Substituting P(x,y, t|x0,y0) as given in Equation (3.18) into this equation,Γ(t|x0) = 14Dpit∫ ∞−∞∫ ∞0[exp(− [(x− x0)2 + (y− y0)2]4Dt)−exp(− [(x+ x0)2 + (y− y0)2]4Dt) ]dxdy .24Evaluating these integrals, we obtainΓ(t|x0) = 12[erfc(− x0√4Dt)− erfc(x0√4Dt)]. (3.19)Thus, this is the probability distribution function for finding the particle in the righthalf-plane at time t, while 1−Γ(t|x0) is the probability that the particle has beenabsorbed by the boundary at time t.(a) Probability of finding the particlein the right half-plane at time t.(b) Probability that the particle has been ab-sorbed by the boundary at time t.Figure 3.9: The probability of finding a particle performing a Brownian mo-tion in the right-half plane at time t given that it started at the point x0,with D = 0.1, and that there is a perfectly absorbing boundary at x = 0.We notice from the plots of the probability distribution function in Equation(3.19) shown in Figure 3.9 that the probability that the diffusing particle is still inthe right half-plane at time t depends on the initial position of the particle. As theparticle starts farther away from the boundary, the time required for the particleto reach the boundary increases, therefore, the probability of finding the particlein the region also increases. We also observe from this figure that the probabilityof finding the particle in the right half-plane at the beginning (t = 0) is 1, and itdecreases with respect to an increase in time. This is because the particle startedits motion in the right half-plane and as time goes on, the chances of it gettingabsorbed by the boundary increases. After a long time, the probability that theparticle is still in the right half-plane tends to zero. This is logical because we25would expect the particle to have been absorbed by the boundary after a long timesince P tends to zero as x tends to infinity. This can easily be seen by taking thelimit of the probability distribution function in Equation (3.19) as t tends to infinity.3.2.2 The disk problemSuppose the particle performing Brownian motion is in a disk-shaped domain of ra-dius a, with a perfectly absorbing boundary at r = a as shown in Figure 3.10. LetΩFigure 3.10: A particle performing Brownian motion in a disk-shaped regionwith a perfectly absorbing boundary.be the interior of the disk-shaped region, that is,Ω= {(r,θ) |0< r ≤ a,0≤ θ ≤ 2pi}.Then the probability that the particle is at position −→x ∈ Ω at time t given that itstarted at point −→x 0 ∈Ω at time t = 0 satisfies∂P∂ t= D∇2P , −→x = (r,θ) ∈Ω, t > 0 ,P(−→x ,0) = δ (−→x −−→x 0) ,P = 0 , on r = a, 0≤ θ ≤ 2pi .(3.20)We shall solve this problem using Green’s function and the method of images26for a disk. First, we rewrite the problem by writing the delta function in the initialcondition as a forcing in the PDE, that is,∂P∂ t=D∇2P+ δ (−→x −−→x 0)δ (t), −→x ∈Ω , t > 0,P = 0, on ∂Ω ,P≡ 0, at t = 0 .(3.21)Let G f be the free space Green’s function for the diffusion equation which satisfies∂G f∂ t=D∇2G f + δ (−→x −−→x 0)δ (t) , −→x ∈Ω , t > 0,G f ≡ 0, at t = 0 .Then from Equation (3.9), we have that the free space Green’s function isG f =14piDtexp(−|−→x −−→x 0|24Dt). (3.22)The idea of the method of images is to find a harmonic function H whichsatisfies the PDE in Equation (3.20), such that P = G f +H, with H = −G f on∂Ω.Let r = |−→x | , r0 = |−→x0 | , r1 = |−→x −−→x0 | , r2 = |−→x0−−→ζ | as shown in Figure 3.11.Then in polar coordinates, Equation (3.22) becomesG f =14piDtexp(− r214Dt). (3.23)We know that r2 = r1(ar)whenever r0 = a. Therefore, we choose the harmonicfunction H to beH=14piDtexp(−Cr224Dt),where C is a constant to be determined.27Figure 3.11: A sketch of the method of images for a disk.Now set P = G f +H and we achieveP =14piDt[exp(− r214Dt)− exp(−Cr224Dt)]. (3.24)Next, we need to find the constant C so that P satisfies the boundary condition.Using the fact that when r0 = a, r2 = r1(ar)and the boundary condition P = 0, weobtained C = r2a2 , and soP =14piDt[exp(− r214Dt)− exp(− r2r224a2Dt)]. (3.25)But from Figure 3.11, using the cosine rule, we haver21 = r2+ r20−2rr0 cos(θ0−θ) ,r22 =1r2(a4+ r2r20−2a2rr0 cos(θ0−θ)) .(3.26)28Therefore,P(r,θ , t) =14piDt[exp(− 14Dt(r2+ r20−2rr0 cos(θ0−θ)))−exp(− 14a2Dt(a4+ r2r20−2a2rr0 cos(θ0−θ)))].(3.27)This is the probability distribution function for finding a particle performing Brow-nian motion in a disk-shaped region of radius a, at a point (r,θ) inside the region,given that it started at a point (r0,θ0) and that there is a perfectly absorbing bound-ary at r = a. We need to integrate this function over the disk-shaped region to getthe probability distribution function of finding the particle inside the region at timet. That is,Γ(t|r0) =∫ΩP(r,θ , t)d−→x =∫ 2pi0∫ a0P(r,θ , t) r dr dθ .Substituting Equation (3.27) into this equation,Γ(t|r0) = 14piDt[∫ 2pi0∫ a0exp(− 14Dt(r2+ r20−2rr0 cos(θ0−θ)))r dr dθ∫ 2pi0∫ a0exp(− 14a2Dt(a4+ r2r20−2a2rr0 cos(θ0−θ)))r dr dθ].Evaluating these integrals using the same approach we used in the case of a disk-shaped region with an imaginary boundary, we haveΓ(t|r0) =(1−Q1(r0√2Dt,a√2Dt))− exp(r20−a24Dt)(1−Q1(ar0√2aDt,a√2aDt)),(3.28)where Q1 is the Marcum Q-function of order one defined in Equation (3.12).Remark 3.2.3. The disk problem for a partially absorbing boundary cannot besolved by the method of images because we cannot find the constant C (in Equa-tion (3.24)) for which P satisfies the Robin boundary condition. Although, othermethods such as separation of variables and the conformal mapping method can29be used to solve the problem.Equation (3.28) is the probability distribution function for finding a particleperforming Brownian motion in a disk-shaped region of radius a, given that itstarted at a distance r0 away from the center of the region and that the regionhas a perfectly absorbing boundary. We notice that the distribution function isindependent of θ . This shows that the probability of finding the particle in theregion is independent of the angle of rotation of the disk. Figure 3.12 shows theplots of the distribution function in Equation (3.28). In this figure, we used a fixedinitial position for the particle and varied the radius of the region. As we wouldexpect, the probability of finding the particle in the disk-shaped region at time t = 0is one since the particle started from inside the region. This probability decreasesover time because the chances of the particle getting absorbed on the boundaryincreases as time goes on.(a) Probability of finding the particleinside the disk-shaped region.(b) Probability that the particle has been ab-sorbed by the boundary at time t.Figure 3.12: The probability of finding a particle performing Brownian mo-tion inside a disk-shaped of radius a at time t, given that it started at adistance r0 = 5 away form the center of the region, with D = 0.1, andthat the region has a perfectly absorbing boundary.We would also expect the probability of finding the particle inside the regionto converge to zero after a long time because the particle is assumed to vanishwhen it hits the boundary (perfectly absorbing boundary), and the chances of theparticle hitting the boundary increases over time. However, the rate at which the30probability converges to zero depends largely on the surface area of the disk-shapedregion. This can be seen in Figure 3.12a.So far, we have derived several probability distribution functions for finding aparticle performing Brownian motion in a region (that is, Equations (3.8), (3.11),(3.19), and (3.28)). These distribution functions can be fit to data to enable one con-clude whether or not an experimental particle of interest is experiencing a restrictedmotion. To do this, we consider the appropriate probability distribution functionsthat best model our experimental set-up and collect data such as the initial positionof the particle x0 (r0 for a disk-shaped region), final position of the particle, itsdiffusion coefficient D, and final experimental time. These pieces of informationcan be used in the probability distribution functions to get the distribution of prob-ability for finding the particle in a region, which can then be used together with thefinal position of the particle to determine whether or not the particle is experiencinga restricted motion.In the next chapter, we derive probability distribution functions for the caseof a partially absorbing boundary and use the distribution function together withmaximum likelihood estimation to estimate the rate of absorption on the boundary.31Chapter 4Partially Absorbing BoundaryIn this chapter, we derive probability distribution functions for finding a parti-cle performing Brownian motion in a bounded domain with a partially absorb-ing boundary. We consider two different geometries: the right half-plane and adisk-shaped region. These distribution functions are used together with maximumlikelihood estimation to estimate the rate of absorption on the boundary. In addi-tion, we use the distribution functions to calculate the mean and variance of firstpassage time for the particle to exit the regions.4.1 Half-plane problemConsider a particle performing Brownian motion in the right half-plane with a par-tially absorbing boundary at x = 0. As illustrated in Figure 4.1, if the diffusingparticle hits the boundary, it is either absorbed by the boundary (in blue) with aspecified probability, or reflected back (in red) into the right half-plane.Erban and Chapman [4] showed that the partially absorbing boundary is accu-rately modelled by a Robin boundary condition. If we let Ω be the right half-plane,that is, Ω = {(x,y)|0≤ x< ∞,−∞< y< ∞}, then the probability of finding theparticle at position−→x ∈Ω at time t, given that it started at position−→x 0 ∈Ω at time32Figure 4.1: A particle performing Brownian motion in the right half-planewith a partially absorbing boundary at x = 0.t = 0, P(−→x , t|−→x 0,0) satisfies∂P∂ t= D∇2P , −→x ∈Ω , t > 0 ,P(−→x ,0) = δ (−→x −−→x 0) ,D∂P∂x= κP , x = 0, −∞< y< ∞,P−→ 0 , as x−→ ∞ ,(4.1)where κ is the rate of absorption on the boundary.Let us rewrite this problem by placing the delta function in the initial conditionas a source function in the PDE,∂P∂ t= D∇2P+δ (x− x0,y− y0)δ (t), 0≤ x< ∞, −∞< y< ∞, t > 0,D∂nP−κP = 0, on x = 0, −∞< y< ∞,P≡ 0, at t = 0.(4.2)33Let P=G f + Ĝ, where G f is the free space Green’s for the diffusion equation, andĜ is a function to be determined. Substituting P=G f + Ĝ into the PDE in problem(4.2), we have(G f + Ĝ)t = D∇2(G f + Ĝ)+δ (x− x0,y− y0)δ (t) .Rearranging,∂G f∂ t−D∇2G f =D∇2Ĝ− ∂ Ĝ∂ t +δ (x− x0,y− y0)δ (t), .Substituting P = G f + Ĝ into the Robin boundary condition in problem (4.2) givesD∂nG f −κG f =−D∂nĜ+κĜ .Therefore, problem (4.2) becomes∂G f∂ t−D∇2G f =D∇2Ĝ− ∂ Ĝ∂ t +δ (x− x0,y− y0)δ (t), 0≤ x< ∞,−∞< y< ∞, t > 0,D∂nG f−κG f =−D∂nĜ+κĜ, on x = 0, −∞< y< ∞,G f ≡ 0, at t = 0 .(4.3)Next, we split this problem into two new problems. The first problem is∂G f∂ t−D∇2G f = δ (x− x0,y− y0)δ (t), −∞< x,y< ∞, t > 0,G f ≡ 0, at t = 0 .(4.4)This is the problem for the free space Green’s function for the diffusion equationand the solution is the Green’s function we derived earlier in Equation (3.7),G f =14piDtexp(− 14Dt[(x− x0)2 + (y− y0)2]). (4.5)34The second problem isD∇2Ĝ− ∂ Ĝ∂ t= 0, 0≤ x< ∞, −∞< y< ∞, t > 0,D∂nG f−κG f =−D∂nĜ+κĜ, on x = 0, −∞< y< ∞, .(4.6)Here, we shall use the idea of the method of images to choose a solution to thePDE in the problem which contains an unknown function called the source densityfunction (see page 478 of [19]). This solution together with the boundary conditionwill be used to find the suitable source density function. We writeĜ =14piDte−14Dt [(x+x0)2+(y−y0)2]+14piDt∫ −x0−∞γ(s)e−14Dt [(x−s)2+(y−y0)2]ds , (4.7)where γ(s) is the source density function, and it is assumed to decay fast enoughat infinity such that the integral converges and differentiation under the integral ispossible. In the Robin boundary condition D∂nG−κG = 0, if κ = 0, we have theNeumann boundary condition, ∂nG = 0, and so we expectĜ =14piDtexp(− 14Dt[(x+ x0)2 + (y− y0)2]).So thatP(x,y, t|x0,y0) = 14piDt exp(− 14Dt[(x− x0)2 + (y− y0)2])+14piDtexp(− 14Dt[(x+ x0)2 + (y− y0)2]),(4.8)which is the probability distribution function for finding the diffusing particle atposition (x,y) at time t when the boundary at x = 0 is perfectly reflecting. There-fore, we want the source density function to be zero when κ = 0. To get this, weneed to find the appropriate γ(s) such that Ĝ satisfies the boundary condition inProblem (4.6).Substituting G f and Ĝ into the boundary condition in problem (4.6) and sim-35plifying, we have− 2κ4piDte−14Dt [x20+(y−y0)2] =κ4piDt∫ −x0−∞γ(s)e−14Dt [s2+(y−y0)2] ds− 14pit∫ −x0−∞γ(s)∂∂x[e−14Dt [(x+s)2+(y−y0)2]]∣∣∣x=0ds .Multiplying through by 4piDt, and using the fact that ∂∂x ≡ − ∂∂ s on the boundary(x = 0), we obtain−2κ e− 14Dt [x20+(y−y0)2] = κ∫ −x0−∞γ(s)e−14Dt [s2+(y−y0)2] ds+D∫ −x0−∞γ(s)∂∂ s[e−14Dt [s2+(y−y0)2]]ds .Using integration by parts for the second integral on the right hand side,−2κ e− 14Dt [x20+(y−y0)2] = κ∫ −x0−∞γ(s)e−14Dt [s2+(y−y0)2] ds+Dγ(−x0)e− 14Dt [x20+(y−y0)2]−D∫ −x0−∞γ ′(s)e−14Dt [s2+(y−y0)2] ds ,(γ ′(s)≡ ∂γ(s)∂ s).Simplifying and collecting the integrals together, we have(2κ +Dγ(−x0))e− 14Dt [x20+(y−y0)2] = (Dγ ′(s)−κ γ(s))∫ −x0−∞e−14Dt [s2+(y−y0)2] ds .We notice that for this to hold, we must have Dγ ′(s)−κ γ(s)= 0 and 2κ +Dγ(−x0)=0, for s<−x0. Therefore, we have the initial value problemDγ ′(s)−κγ(s) = 0 , s<−x0 ,Dγ(−x0) =−2κ .Solving this initial value problem, we obtain γ(s) = −2κD exp( κD(x0+ s)). Substi-tuting for γ(s) in Equation (4.7) givesĜ =14piDte−14Dσ [(x+x0)2+(y−y0)2]− 2κ4piD2t∫ −x0−∞eκD (x0+s) e−14Dt [(x−s)2+(y−y0)2] ds .36Substituting G f and Ĝ into P = G f + Ĝ, we haveP(−→x , t) = 14piDt[e−14Dt [(x−x0)2+(y−y0)2]+ e−14Dt [(x+x0)2+(y−y0)2]]− κ2piD2t∫ −x0−∞eκD (x0+s) e−14Dt [(x−s)2+(y−y0)2] ds .Evaluating the integral, we getP(−→x , t) = 14piDt[e−14Dt [(x−x0)2+(y−y0)2]+ e−14Dt [(x+x0)2+(y−y0)2]]− κpiD(x+ x0)+2piDtκ(e−14Dt [(x+x0)2+(y−y0)2]).(4.9)This is the probability distribution function for finding a particle performing Brow-nian motion in the right half-plane at position−→x at time t, given that it started fromposition −→x 0 at time t = 0 and that there is a partially absorbing boundary at x = 0.We shall integrate this function over the right half-plane to obtain the probabil-ity distribution function for finding the particle in the right half-plane at any giventime t, that isΓ(t|x0) =∫ ∞−∞∫ ∞0P(x,y, t|x0,y0) dxdy .Substituting P(−→x , t) as given in Equation (4.9) into this equation,Γ(t|x0) = 14piDt∫ ∞−∞∫ ∞0[e−14Dt [(x−x0)2+(y−y0)2]+ e−14Dt [(x+x0)2+(y−y0)2]]dxdy−∫ ∞−∞∫ ∞0κpiD(x+ x0)+2piDtκ[e−14Dt [(x+x0)2+(y−y0)2]]dxdy .Evaluating the integrals in the first term, we haveΓ(t|x0) = 12[erfc(− x0√4Dt)+ erfc(x0√4Dt)]−∫ ∞0√4piDt κpiD(x+ x0)+2piDtκ(e−14Dt (x+x0)2)dx .37Since 12[erfc(− x0√4Dt)+ erfc(x0√4Dt)]= 1, we can writeΓ(t|x0) = 1−∫ ∞0√4piDt κpiD(x+ x0)+2piDtκ(e−14Dt (x+x0)2)dx .Let Γ(t|x0) = 1−√4piDt κ I, whereI =∫ ∞01piD(x+ x0)+2piDtκ(e−14Dt (x+x0)2)dx .Rewriting I, we haveI =1piD∫ ∞0(1x+(x0+2tκ))e−14Dt (x+x0)2dx .We notice that this integral does not converge, therefore, to evaluate the integralanalytically, we Taylor expand the fraction in the integrand near x = 0 to haveI ' 1Dpi∫ ∞0(1(x0+2tκ)− x(x0+2tκ)2+x2(x0+2tκ)3+ . . .)e−14Dt (x+x0)2dx .I ' 1Dpi(x0+2tκ)[∫ ∞0e−14Dt (x+x0)2dx− 1(x0+2tκ)∫ ∞0xe−14Dt (x+x0)2dx+ . . .].Evaluating the integrals, we obtainI ' 1Dpi(x0+2tκ)[√4piDt2erfc(x0√4Dt)− 12(x0+2tκ)(4Dt e(−x204Dt )−√4piDt x0 erfc(x0√4Dt))+ . . .].Therefore,I ' 1Dpi(x0+2tκ)[√4piDt2erfc(x0√4Dt)− 12(x0+2tκ)(4Dt e(−x204Dt )−√4piDt x0 erfc(x0√4Dt)) ].(4.10)38It is important to point out that this solution is valid for 1< |x0+2κt|. ThusΓ(t|x0) = 1−√4piDt κpiD(x0+2tκ)[√4piDt2erfc(x0√4Dt)− 12(x0+2tκ)(4Dt exp(− x204Dt))−√4piDt x0 erfc(x0√4Dt))].Simplifying, we have the probability distribution function for finding the particlein the right half-plane at time t asΓ(t|x0) = 1−[(1+x0(x0+2tκ))2t κ(x0+2tκ)erfc(x0√4Dt)− 4(Dt)3/2κ√piD(x0+2tκ)2exp(− x204Dt)].(4.11)where κ is the rate of absorption on the boundary, x0 is the initial position of theparticle, and D is its diffusion coefficient.We observe from the Robin boundary condition D ∂P∂x = κP that if κ = 0, theboundary condition becomes a perfectly reflecting (Neumann) boundary condition,∂P∂x = 0. And if we set κ = 0 in the probability distribution function in Equation(4.11), we have Γ(t|x0) = 1 which corresponds to the probability that the particleis in the right half-plane at time t, given that the boundary at x = 0 is perfectlyreflecting. This is logical because the particle will remain in the region for all timesince the boundary is perfectly reflecting.On the other hand, if we take the limit of the Robin boundary condition as κtends to infinity, the boundary condition becomes the perfectly absorbing boundarycondition, P = 0. Also, taking the limit of the distribution function in Equation(4.11) as κ tends to infinity, we haveΓ(t|x0) = 1− erfc(x0√4Dt)=12[erfc(− x0√4Dt)− erfc(x0√4Dt)],which is the same as the probability distribution function derived in Chapter 3(Equation (3.19)) for finding a diffusing particle in the right half-plane, where theboundary at x = 0 is perfectly absorbing. This shows that the probability distri-39bution function for a partially absorbing boundary is a weighted average of thedistribution functions for an absorbing boundary and that of a reflecting boundary.Next, we simulate trajectories with the particle based simulation software [2],calculate the probability of finding a diffusing particle in the right half-plane fromthe simulation, and compare the result to that of the probability distribution func-tion in Equation (4.11). To calculate the probability of finding a particle performingBrownian motion in a region at a specific time with simulation, we simulate severaltrajectories that all start at the same point, count the number of trajectories that areremaining in the region at each time step and divide by the initial number of tra-jectories. This gives the probability that a single trajectory that started a Brownianmotion at that point is still in the region at each time step. We specify a particularrate of absorption on the boundary κ , and use the relation in Equation (2.9) to getthe corresponding probability of absorption on the boundary which is then used forthe simulation.(a) Probability of finding the particlein the right half-plane.(b) Absolute difference of the analytic andsimulated result.Figure 4.2: The comparison of the probability of finding the particle in theright half-plane obtained from the simulation and that of the probabilitydistribution function in Equation (4.11).Figure 4.2a shows a comparison of the probability that the particle perform-ing Brownian motion in the right half-plane is still in the region at time t, calcu-lated from simulation, and that of the probability distribution function in Equation(4.11), while Figure 4.2b shows the absolute difference of the two results. For the40result obtained from simulation, we averaged over 10 different simulations using10,000 particles for each simulation, with κ = 0.5 which corresponds to probabil-ity of absorption P= 0.044, D = 0.1, and x0 = 2. We can see from this figure thatthe two result agree to a large extent.Let us calculate the mean and variance of first passage time for the diffusingparticle to exit the right half-plane. First, we notice that the probability that theparticle is still in the right half-plane at time t given that it started at point x0 isequivalent to the probability that the time it takes the particle to be absorbed by theboundary given that it started from x0 is greater than t, that is,Γ(t|x0) =∫ ∞tΨ(τ|x0)dτ , (4.12)where Ψ(τ|x0) is the probability distribution function for the first passage time forthe particle to exit the region. From Equation (4.12), we haveΨ(t|x0) =− ∂∂ tΓ(t|x0) . (4.13)Substituting Γ(t|x0) as given in Equation (4.11) into Equation (4.13), we obtainΨ(t|x0) = 4κx20(x0+2tκ)3erfc(x0√4Dt)+[x20+2tκ(x0− 3Dκ +8Dt(x0+2tκ))]× κ√Dpit(x0+2tκ)2exp(− x204Dt).(4.14)Let Π(t|x0) be the mean first passage time for the particle to exit the region, giventhat it started at x0. Since the first moment of the probability distribution for thefirst passage time is the mean first passage time, we haveΠ(t|x0) =∫ ∞0tΨ(t|x0)dt . (4.15)The second moment of the distribution isΦ(t|x0) =∫ ∞0t2Ψ(t|x0)dt . (4.16)41Therefore, the variance of first passage time for the particle to exit the region giventhat it started at x0 isV(t|x0) =Φ(t|x0)−Π(t|x0)2 . (4.17)We substitute Equation (4.14) into Equations (4.15) and (4.16), and evaluate theresulting integrals numerically. These results are then used to compute the meanand variance of first passage time for the particle to exit the right half-plane.(a) Mean first passage time. (b) Variance of first passage time.Figure 4.3: The mean and variance of first passage time for a particle per-forming Brownian motion in the right half-plane, with a partially ab-sorbing boundary at x = 0.Figure 4.5 shows the plots of the mean and variance of first passage time for adiffusing particle to exit the right half-plane for different values of κ . We observefrom Figure 4.5a that the mean first passage time increases as the initial position ofthe particle moves farther away from the boundary. This is not surprising becausethe farther away from the boundary the particle starts, the longer it takes it to reachthe boundary. Also we observe that as the rate of absorption on the boundary in-creases, the mean first passage time decreases. This is because the particle is easilyabsorbed on the boundary as the rate of absorption on the boundary increases. Themean and variance of first passage time tells us how long we should expect a parti-cle performing Brownian motion in the right half-plane to stay in the region beforeit gets absorbed by the boundary.424.2 The disk problemSuppose the diffusing particle is in a disk-shaped region of radius a with a partiallyabsorbing boundary. Let Ω= {(r,θ) |0< r ≤ a,0≤ θ ≤ 2pi}, then the probabilityof finding the particle at position−→x =(r,θ) at time t, given that it started at position−→x 0 = (r0,θ0) at time t = 0 satisfies∂P∂ t= D∇2P , −→x = (r,θ) ∈Ω , t > 0 ,D∂P∂n=−κP , on ∂Ω ,P is finite as r −→ 0 ,P(−→x ,0) = δ (−→x −−→x 0) .(4.18)where κ is the rate of absorption on the boundary, and D is the diffusion coefficientof the particle.We shall solve this problem using separation of variables. First, let us rewritethe problem by writing the delta function in the initial condition as a forcing in thePDE, that is,∂P∂ t=D(1r∂P∂ r+∂ 2P∂ r2)+δ (r− r0)δ (t)2pir, 0< r ≤ a , t > 0 ,D∂P∂ r=−κP , on r = a ,P is finite as r −→ 0 ,P≡ 0, at t = 0 .(4.19)It is important to note that we have assumed that there is uniformity in the θdirection so that the problem is independent of θ .From the problem in Equation (4.19), we consider the homogeneous part of thePDE∂P∂ t=D(1r∂P∂ r+∂ 2P∂ r2), (4.20)43together with the following boundary conditionsD∂P∂ r=−κP , on r = a ,P is finite as r −→ 0 .(4.21)Let P(r, t) = R(r)T (t), then from Equations (4.20) and (4.21), we haveR(r)ddtT (t) = D(1rT (t)ddrR(r)+T (t)d2dr2R(r)),DddrR(r) =−κR(r) on r = a ,R(r) is finite as r −→ 0 .(4.22)Let us consider1DT (t)ddtT (t) =1R(r)(1rddrR(r)+d2dr2R(r)).For this to hold, it must equal a constant, that is1DT (t)ddtT (t) =1R(r)(1rddrR(r)+d2dr2R(r))=−λ 2 . (4.23)where λ is a constant.Therefore, we have the following equations1DT (t)ddtT (t) =−λ 2 and 1R(r)(1rddrR(r)+d2dr2R(r))=−λ 2 . (4.24)Consider1R(r)(1rddrR(r)+d2dr2R(r))=−λ 2 .Multiplying through by r2R(r) and re-arranging, we have the Bessel equationr2R′′(r)+ rR′(r)+ r2λ 2R(r) = 0 .44And the solution of this equation isR(r) = AJ0(λ r)+BY0(λ r) ,where J0(z) is the Bessel function of the first kind of order zero, and Y0(z) is theBessel function of the second kind of order zero.To satisfy the condition that R(r) is finite as r tends to zero, we must have B= 0because Y0(λ r) is unbounded as r tends to zero. Therefore,R(r) = AJ0(λ r) .Imposing the second boundary condition in Equation (4.22), we haveDddrJ0(λ r) =−κJ0(λ r) , on r = a .And this gives the following equation for the eigenvalues of the systemκJ0(λa)+DλJ1(λa) = 0. (4.25)Let βm be the mth zero of the eigenvalue equation, then λm = βma . Therefore, wehave the eigen conditionaκJ0(βm)+DβmJ1(βm) = 0, m = 1,2,3, . . . (4.26)Without loss of generality, let Am = 1 for all m. ThenRm(r) = J0(λmr) = J0(βmar), m = 1,2,3, . . . (4.27)Recall that P(r, t) = R(r)T (t). This implies that Pm(r, t) = Rm(r)Tm(t) for all m.Therefore, using the principle of superposition, we haveP(r, t) =∞∑m=1Tm(t)Rm(r) =∞∑m=1Tm(t)J0(βmar). (4.28)45Now, let∞∑m=1Hm(t)Rm(r) =δ (r− r0)δ (t)2pir, (4.29)Substituting Rm(r) into this equation, and multiplying through by r, we have∞∑m=1r Hm(t)J0(βmar)=δ (r− r0)δ (t)2pi.Multiplying both sides by J0(βna r), and integrating with respect to r from 0 to a,∞∑m=1Hm(t)∫ a0r J0(βmar)J0(βnar)dr =δ (t)2pi∫ a0J0(βnar)δ (r− r0)dr .By the orthogonality property of Bessel functions,∫ a0r J0(βmar)J0(βnar)dr = 0 , m 6= n .And so, we haveHm(t)∫ a0r J20(βmar)dr =δ (t)2piJ0(βmar0).But ∫ a0r J20(βmar)dr =a22[J20 (βm)+ J21 (βm)].Therefore,Hm(t) =J0(βma r0)pia2[J20 (βm)+ J21 (βm)]δ (t) . (4.30)From Equation (4.19), we have∂P∂ t=D(1r∂G∂ r+∂ 2G∂ r2)+δ (r− r0)δ (t)2pir, (4.31)46Substituting Equations (4.28) and (4.29) into Equation (4.31),∞∑m=1T ′m(t)Rm(r) = D[1r∞∑m=1Tm(t)R′m(r)+∞∑m=1Tm(t)R′′m(r)]+∞∑m=1Hm(t)Rm(r) ,= D[∞∑m=1Tm(t)(1rR′m(r)+ R′′m(r))]+∞∑m=1Hm(t)Rm(r) .(4.32)From the Bessel equation, we haveR′′m(r)+1rR′m(r) =−λ 2mRm(r) .Therefore, Equation (4.32) becomes∞∑m=1T ′m(t)Rm(r) =∞∑m=1(−Dλ 2mTm(t)+Hm(t))Rm(r) .And from this equation, we haveT ′m(t) =−Dλ 2mTm(t)+Hm(t) .Substituting Hm(t) as given in Equation (4.30) into this equation,T ′m(t)+Dλ2mTm(t) =J0(βma r0)pia2[J20 (βm)+ J21 (βm)]δ (t) .This is a first order linear ODE. We shall solve this ODE using the integratingfactor method. The integrating factor is exp(Dλ 2mt), therefore,ddt[exp(Dλ 2mt)Tm(t)]=J0(βma r0)pia2[J20 (βm)+ J21 (βm)]δ (t)exp(Dλ 2mt) .47Integrating both sides of the equation with respect to t,exp(Dλ 2mt)Tm(t) =J0(βma r0)pia2[J20 (βm)+ J21 (βm)] ∫ t0δ (s)exp(Dλ 2ms)ds+Em , (4.33)where Em are constants of integration.We shall use the initial condition P(r,0) = 0 to find the constants Em. We haveP(r,0) =∞∑m=1Tm(0)Rm(r) .Setting t = 0 in Equation (4.33), we get Tm(0) = Em. And substituting this intoP(r,0) givesP(r,0) =∞∑m=1Tm(0)Rm(r) =∞∑m=1Em Rm(r) .Using the orthogonality property of Bessel functions and the fact that P(r,0) = 0,we have that Em = 0 for all m. Therefore, from Equation (4.33) we obtainTm(t) =J0(βma r0)pia2[J20 (βm)+ J21 (βm)] exp(−Dλ 2mt) .Substituting λm = βma and J1(βm) =κDλmJ0(βm), we haveTm(t) =D2β 2ma2piJ0(βma r0)(D2β 2m+κ2a2)J20 (βm)exp(−Dβ2ma2t). (4.34)Substituting Tm(t) and Rm(r) into Equation (4.28), we haveP(r, t|r0) = D2a2pi∞∑m=1β 2mJ0(βma r0)J0(βma r)(D2β 2m+κ2a2)J20 (βm)exp(−Dβ2ma2t). (4.35)This is the probability distribution function for finding a particle performing Brow-nian motion in a disk-shaped region of radius a, at a distance r from the center ofthe region at time t, given that it started at a distance r0 from the center, and that48the region has a partially absorbing boundary.Lastly, we integrate Equation (4.35) to get the probability distribution functionfor finding the particle inside the disk-shaped region at time t.Γ(t|r0) = D2a2pi∞∑m=1β 2m J0(βma r0)(D2β 2m+κ2a2)J20 (βm)exp(−Dβ2ma2t)∫ 2pi0∫ a0r J0(βmar)dr dθ .Γ(t|r0) = 2D2∞∑m=1βm J0(βma r0)J1(βm)(D2β 2m+κ2a2)J20 (βm)exp(−Dβ2ma2t). (4.36)Figure 4.4a shows the comparison of the probability of finding a diffusing par-ticle in a disk-shaped region of radius a = 5 calculated from simulation and that ofthe probability distribution function in Equation (4.36), while Figure 4.4b showsthe absolute difference of the two results. For the simulated result, we averagedover 10 different simulations each with 10,000 particles starting from a distancer0 = 2 away the center of the region, with diffusion coefficient D = 0.1, and therate of absorption on the boundary κ = 0.85. We observe from theses figures thatthe simulated result is in accordance with that obtained from the probability distri-bution function.(a) Probability of finding the particle inthe disk-shaped region.(b) Absolute difference of the analytic andnumerical result.Figure 4.4: The comparison of the probability of finding a particle in the disk-shaped region calculated from simulation and that of the probabilitydistribution function in Equation (4.36).49Similar to the right half-plane problem, we want to calculate the mean andvariance of first passage time for a particle performing Brownian motion in a disk-shaped region to exit the region. Substituting Equation (4.36) into Equation (4.13),we obtain the distribution function for first passage time for a diffusing particle inthe disk-shaped region,Ψ(t|r0) = 2D3a2∞∑m=1β 3m J0(βma r0)J1(βm)(D2β 2m+κ2a2)J20 (βm)exp(−Dβ2ma2t). (4.37)Using Equations (4.15), (4.16), and (4.17), we can get the mean and variance offirst passage time for the particle. We substitute Equation (4.37) into the integrals inEquations (4.15) and (4.16), and evaluate the resulting integrals numerically to getthe mean and the second moment of the distribution. These results are then usedin Equation (4.17) to compute the variance of first passage time for the particle.Figure 4.5 shows the plots of the mean and variance of first passage time for aparticle performing Brownian motion in a disk-shaped region of radius a = 7, witha partially absorbing boundary.(a) Mean first passage time. (b) Variance of first passage time.Figure 4.5: The mean and variance of first passage time for a particle per-forming Brownian motion in a disk-shaped region with a partially ab-sorbing boundary.For this figure, we used D= 0.1, and varied the rate of absorption on the bound-ary and the initial position of the particle. It is important to note that the initial50position of the particle r0 is the distance of the particle from the centre of the disk-shaped region. As we would expect, the mean first passage time for the particleto exit the region decreases as the initial position of the particle gets closer to theboundary. This is due to the fact that the time it takes the particle to reach theboundary when it starts close to the boundary is less compared to when the particlestarts farther away from the boundary. Also, like in the right half-plane problem,we observe that the mean first passage time decreases as the rate of absorption onthe boundary increases. This is because the particle is easily absorbed as the rateof absorption on the boundary gets larger. Unlike the case of a particle in the righthalf-plane where the mean and variance of first passage time are linear with respectto the initial position of the particle, they are non-linear for this case, and this is asa result of the geometry of the region.4.2.1 Estimating the rate of absorption on the boundaryIn this section, we present a mathematical technique for estimating the rate of ab-sorption on a partially absorbing boundary using the probability distribution func-tions we derived for a particle performing Brownian motion in regions with par-tially absorbing boundary together with maximum likelihood estimation.First, let us derive the likelihood function. We define the likelihood function ofhaving the parameter κ given some observations as`(κ|Observations) = P(Observations;κ) . (4.38)where P(Observations;κ) is the probability of having these observations with pa-rameter κ . There are only two possible observations here;• the diffusing particle has been absorbed by the boundary at time t• the particle is still in the region at time tLet Σ be our region of interest. Suppose there are n+m identical particlesstarting a Brownian motion at a specific point in Σ at time t = 0, the motion ofeach particle is independent. Let P(Particle in Σ;κ) be the probability that a par-ticle is still in the region of interest at time t given that κ is a parameter, and letP(Particle absorbed;κ) be the probability that a particle has been absorbed by the51boundary at time t. Then, the probability of having the observations that n parti-cles are still remaining in the region, while m of them have been absorbed by theboundary at time t, given that κ is a parameter isP(Observations;κ) = P(Particle is in Σ;κ)n×P(Particle absorbed;κ)m .But P(Particle absorbed;κ) = 1−P(Particle is in Σ;κ), therefore,P(Observations;κ) = P(Particle is in Σ;κ)n× (1−P(Particle is in Σ;κ))m .And so, the likelihood function is`(κ|Observations) = P(Particle is in Σ;κ)n× (1−P(Particle is in Σ;κ))m ,where n is the number of particles left in the region of interest, and m is the numberof particles that were absorbed by the boundary after a specific time.Thus, the likelihood function can be written as`(κ|Observations) = (Γ(t|x0))n× (1− (Γ(t|x0))m , (4.39)where Γ(t|x0) is the probability distribution function for finding a particle in theregion of interest at time t given that it started from the point x0 in the region. Thatis, Equation (4.11) for the right half-plane problem, and Equation (4.36) for thedisk problem.We need to maximize the likelihood function with respect to κ , and the κvalue that maximizes this function is the most probable rate of absorption on theboundary with respect to the observations or data.To validate this idea, we simulate trajectories using the particle based simula-tion software [2] with a specific rate of absorption on the boundary, collect datafrom the simulation in terms of the number of trajectories remaining in the regionof interest and those that have been absorbed by the boundary after a specifiedtime. The data collected are used in the likelihood function in Equation (4.39), andthen the resulting function is maximized with respect to κ to recover the rate ofabsorption used for the simulation.52(a) (b)Figure 4.6: Estimates of the rate of absorption on the boundary with stan-dard deviation error bars based on simulated data considering differentparticle initial positions x0, for the half-plane problem.Figure 4.6 shows the result of our estimates of the rate of absorption on theboundary with standard deviation error bars, based on simulated data for a particleperforming Brownian motion in the right half-plane. Here, we used a fixed rateof absorption κ = 0.75, diffusion coefficient D = 0.1, and vary the initial positionof the particles. For each initial position, we simulate 10,000 trajectories for timet = 20, and then use our estimation technique to recover the rate of absorption onthe boundary using the simulated data. The estimates in Figure 4.6a were obtainedwith 20 simulations, while those in Figure 4.6b were obtained with 30 simulations.We observe from these figures that we have good estimates, and that the accuracy ofour estimate increases as the number of simulations used increases. This suggeststhat more identical experimental procedures should be considered when estimatingthe rate of absorption on the boundary in order to increase the accuracy of theestimate. In addition, we notice that the accuracy of the estimates decreases as theinitial position of the particles moves farther away from the boundary. This is dueto the fact that when the particles start far away from the boundary, it takes moretime for them to get to the boundary and they spread out more compared to whenthe start closer to the boundary. To get a better estimate when the particles startsfar away from the boundary, we can either increase the final time for the simulation53or increase the number of particle in the simulation.(a) (b)Figure 4.7: Estimates of the rate of absorption on the boundary with standarddeviation error bars based on simulated data using different κ values,for the half-plane problem.In Figure 4.7, we show our estimates of the rate of absorption on the boundarywith standard deviation error bars for a situation where we fixed the initial posi-tion of the particle and diffusion coefficient, but vary the rate of absorption on theboundary. The rates of absorption used are 0.5,0.8,1.1, and 1.5, with diffusion co-efficient D= 0.1, final time t = 20, and particle initial position x0 = 2. The result inFigure 4.7a was obtained using 20 different simulations, while that of Figure 4.7bwas obtained using 30 simulations. Similar to the estimates in Figure 4.6, we ob-serve that the estimates get better as the number of simulations increases, althoughthis comes at the expense of computational time.Next, we present similar results for the disk problem. Figure 4.8 shows theestimates of the rate of absorption on the boundary of the disk-shaped region withstandard deviation error bar. For these estimates, we consider a region of radiusa = 5, with rate of absorption on the boundary κ = 0.5, diffusion coefficient D =0.1, final time t = 20, and vary the initial position of the particles. We used 10different simulations to obtain the estimates in Figure 4.8a and 20 simulations forthose in Figure 4.8b. We notice from this figure that we also have better estimateswhen the particle starts close to the boundary and when we use more simulations54like in the case of the half-plane problem.(a) (b)Figure 4.8: Estimates of the rate of absorption on the boundary with stan-dard deviation error bars based on simulated data considering differentparticle initial positions r0, for the disk problem.(a) (b)Figure 4.9: Estimates of the rate of absorption on the boundary with standarddeviation error bars based on simulated data using different κ values,for the disk problem.Figure 4.9a shows the estimates of the rate of absorption on the boundary ofa disk-shaped region with standard deviation error bar for a case where we variedthe rate of absorption on the boundary. We used the same value of the parameters55a, D, and t, as that of Figure 4.8, with initial position r0 = 2, and vary the rate ofabsorption on the boundary. The κ values used for this figure are 0.5,0.8,1.1, and1.5. Similar to the right half-plane problem, we observe that the accuracy of ourestimates decreases as the rate of absorption on the boundary increases.In the next chapter, we shall consider a situation where the rate of absorptionon the boundary fluctuates over time.56Chapter 5Fluctuating BoundaryIn this chapter, we consider scenarios where the boundary fluctuates from oneboundary type to another. This chapter is motivated by a situation where the barriercausing the boundary is made up of some other diffusing particles. We present atechnique for estimating the rate of absorption on the boundary by approximatingthe fluctuating boundary with a specific partially absorbing boundary. Further-more, we consider a situation where the boundary switches between two boundarytypes at a specified time interval, and present mathematical relations for finding theeffective rate of absorption on the boundary as a function of the switching time.5.1 Estimating the rate of absorption on a fluctuatingboundarySo far we have been discussing boundaries without taking a close look at the ma-terial that forms the boundaries. We know that in some cases, the barrier causingthe boundary is not fixed. For example in cellular biology, molecules diffusing ata slow rate, and with high density may restrict the motion of some other moleculesdiffusing at a higher rate, thereby forming a boundary for these molecules. Figure5.1a shows an illustration of a particle performing Brownian motion in the righthalf-plane (in black), where some other particles also performing Brownian mo-tion (in blue) create a barrier forming a boundary at x = 0, and Figure 5.1b showsan illustration of a particle performing Brownian motion in a disk-shaped region57bounded by some particles that are also diffusing.(a) The right half-plane. (b) A disk-shaped region.Figure 5.1: A particle performing Brownian motion in some region boundedby some other diffusing particles.In this situation, we do not expect the boundary to be perfectly absorbing (ortransmitting) or perfectly reflecting for all time, rather we would expect it to fluctu-ate between perfectly reflecting and partially absorbing over time. This boundarycan also be seen as a stochastic switching boundary (see [8]) where the switchingrate depends on the diffusion coefficient of the particles that create the barrier. Animportant question to ask is: How can we estimate the rate of absorption on theboundary?Since the boundary fluctuates between partially absorbing and perfectly reflect-ing, we can model this boundary with a time varying Robin boundary condition,that is, we let the rate of absorption in the Robin boundary condition depends ontime, as shown below,D∂P∂x= κ(t)P , (5.1)where κ(t) is the time-dependent rate of absorption on the boundary.Note that we are only considering a situation where the particle of interest(black in Figure 5.1) is either absorbed by the boundary or reflected back into thedomain. Also, we assume that the boundary can never be free of the molecules thatcreate the barrier, that is , the boundary is never perfectly absorbing.58The challenge with this problem is solving the PDE in Equation (4.1) togetherwith the time varying Robin boundary condition analytically. As a result of this,we consider a much simpler Robin boundary condition using mean field approxi-mation,D∂P∂x= κP , (5.2)where κ is the average of the rates of absorption on the boundary over time.Now we can approximate the fluctuating boundary with a partially absorbingboundary where the rate of absorption on the boundary is the average of the fluc-tuating rates of absorption on the boundary.For the half-plane problem, we have Ω = {(x,y) |0≤ x< ∞, −∞< y< ∞}.The probability that the particle of interest is at position −→x at time t given that itstarted at −→x 0 at time t = 0, P(−→x , t|−→x 0,0) satisfies the following problem;∂P∂ t= D∇2P , −→x ∈Ω , t > 0 ,P(−→x ,0) = δ (−→x −−→x 0) ,D∂P∂x= κP , on x = 0, −∞< y< ∞,P−→ 0 , as x−→ ∞ .(5.3)Notice that this problem is exactly the same with the problem in Equation (4.1)if we replace κ with κ . Therefore, we adopt the solution in Equation (4.11), andso we have that the probability distribution function for finding the particle in theright half-plane at time t given that it started at x0 isΓ(t|x0) = 1−[(1+x0(x0+2tκ))2t κ(x0+2tκ)erfc(x0√4Dt)− 4(Dt)3/2κ√piD(x0+2tκ)2exp(− x204Dt)],(5.4)where κ is the average of the rates of absorption on the boundary over time.Similarly, for a disk-shaped region, Ω = {(r,θ) |0< r ≤ a,0≤ θ ≤ 2pi}, andthe probability of finding the particle of interest at position −→x = (r,θ) at time t59given that it started at position −→x 0 = (r0,θ0) at time t = 0 satisfies∂P∂ t= D∇2P , −→x = (r,θ) ∈Ω , t > 0 ,D∂P∂n=−κP , on ∂Ω ,P is finite as r −→ 0 ,P(−→x ,0) = δ (−→x −−→x 0) .(5.5)From Equation (4.36), the probability distribution function for finding the particlein the disk-shaped region at t given that it started at a distance r0 from the center ofthe region isΓ(t|r0) = 2D2∞∑m=1βm J0(βma r0)J1(βm)(D2β 2m+κ2a2)J20 (βm)exp(−Dβ2ma2t). (5.6)where κ is the average of the rates of absorption on the boundary over time.Next, we verify that the mean field approximation works for these problemsby simulating trajectories in a region where the rate of absorption on the bound-ary fluctuates over time, and then calculate the probability of finding a particleperforming Brownian motion inside the region from the simulation. In our simula-tion, we let the switching on the boundary happen at every unit time interval, andthen calculate the average rate of absorption on the boundary which will be used inthe probability distribution functions in Equations (5.4) and (5.6).Figure 5.2 shows the fluctuation in the rate of absorption on the boundary overtime, and the comparison of the probability for finding the particle of interest inthe right half-plane calculated from simulation and the result obtained using meanfield approximation, while Figure 5.3 shows similar results for a disk-shaped re-gion of radius a= 5. For these figures, we averaged over 10 simulations, each with10,000 trajectories, with diffusion coefficient D = 0.1, x0 = 5 for the half-planeproblem, and r0 = 2 for the disk problem. Although the rates of absorption on theboundaries of the two regions fluctuate over time in the simulations, the probabil-ity of finding the particle of interest in the regions obtained using the mean fieldapproximation still agrees with the probability calculated form simulation. This60(a) Fluctuation on the boundary in termsof κ .(b) Prob. of finding the particle in the righthalf-plane.Figure 5.2: Comparison of the probability of finding the particle of interestin the right half-plane calculated from the simulation and the predictionfrom the mean field approximation.(a) Fluctuation on the boundary in termsof κ .(b) Prob. of finding the particle in the disk-shaped region.Figure 5.3: Comparison of the probability of finding the particle of interestin the a disk-shaped region calculated from the simulation and the pre-diction from the mean field approximation.shows that the mean approximation works for these problems and that the effectiverate of absorption on the boundary is the average of the fluctuating rates of absorp-tion. It is also true that this mean field approximation works in a situation wherethe switching on the boundary happens at random.615.2 Switching boundaryIn this section, we consider situations where the boundary switches between twospecific boundary types. First, we consider a case where the boundary switches be-tween a perfectly reflecting boundary and a specific partially absorbing. Followedby the case where the boundary switches between a perfectly absorbing boundaryand a perfectly reflecting boundary. For both cases, we approximate the switch-ing boundary with a partially absorbing boundary, where the rate of absorption onthe boundary is the effective rate of absorption that best approximates the solution,given the switching on the boundary. We also derive an empirical mathematicalrelation for the effective rate of absorption as a function of the boundary switchingtime.5.2.1 Switching between perfectly reflecting and partially absorbingSuppose the boundary of the region where the particle is performing Brownianmotion switches between a perfectly reflecting boundary and a specific partiallyabsorbing boundary at some specified time interval. Since a perfectly reflectingboundary condition can be obtained from the Robin boundary condition by settingthe rate of absorption on the boundary to be zero (i.e κ = 0), we approach thisproblem from the point of view of a partially absorbing boundary with rate of ab-sorption that switching between two specific rates (say, κ1 = 0 and κ2). We beginby setting the effective rate of absorption on the boundary to be the average of thetwo rates of absorption, that is, κe f f = κ . To validate this approximation, we sim-ulate trajectories in a region where the rate of absorption on the boundary switchesbetween these two rates, and from the simulation, we calculate the probability offinding the particle in the region at each time step and compare with the result ofthe appropriate probability distribution function.Since we are approximating the switching boundary with a partially absorbingboundary whose rate of absorption is κe f f , the resulting problem is the exactly thesame as problem (4.1) for the half-plane problem and problem (4.18) for the diskproblem with κ replaced by κe f f . Therefore, the probability distribution function62for finding the particle in the right half-plane is,Γ(t|x0) = 1−[(1+x0(x0+2tκe f f ))2t κe f f(x0+2tκe f f )erfc(x0√4Dt)− 4(Dt)3/2κe f f√piD(x0+2tκe f f )2exp(− x204Dt)].(5.7)While that of a disk-shaped region is,Γ(t|r0) = 2D2∞∑m=1βm J0(βma r0)J1(βm)(D2β 2m+κ2e f f a2)J20 (βm)exp(−Dβ2ma2t). (5.8)Figure 5.4 shows the comparison of the probabilities calculated from simula-tion and that of the probability distribution function in Equation (5.7).We notice from these plots that irrespective of the value of κ2, there is alwaysa discrepancy between the result from the probability distribution function and thatof the simulation. This suggests that the average of the rates of absorption on theboundary is not the effective rate of absorption. We have presented the plots forthe right half-plane problem only. The analogous results for the disk problem havea similar discrepancy.The next question is; How can we find the effective rate of absorption on theboundary? To answer this question, we fit the probability distribution functionin Equation (5.7) to the probability of finding the particle in the right half-planecalculated from simulation as shown is Figure 5.5. The κ value that best fits thedistribution function to the result from simulation is the effective rate of absorptionon the boundary κe f f . Following the same procedure, we can fit the probabilitydistribution function in Equation (5.8) to simulated results to get the effective rateof absorption on the boundary for the case of a disk-shaped region.Our goal is to derive a mathematical relation that predicts the effective rate ofabsorption κe f f as a function of the rate of absorption κ2 and the switching time ofthe boundary. We define the switching time as a regular time interval at which theboundary switches from one boundary type to the other, and denote it by α . Sincewe can get the effective rate of absorption on the boundary by fitting a probabil-ity distribution function to the probability for finding the particle in the region of63Figure 5.4: Comparison of the probability of finding the particle of interestin the right half-plane calculated from simulation and that of the proba-bility distribution function in Equation (5.7).interest calculated from simulation, we begin by doing this for different values ofα . The values of α considered are between 1/2 and 4 with an increment of 1/2.Figure 5.6 shows the estimates of the effective rate of absorption obtained by fittingthe probability distribution functions in Equations (5.7) and (5.8) to the probabilitycalculated from simulation for different rate of absorption κ2 and switching time.For the estimates in this figure, we averaged over 10 different simulations eachwith 10,000 trajectories, diffusion coefficient D = 0.1, initial position x0 = 5 forthe right half-plane problem, and r0 = 2 for the disk problem with a disk-shaped re-gion of radius a= 5. We notice from this figure that the estimates for the half-planeproblem is similar to that of the disk problem. Also, we observe that for each valueof α , the effective rate of absorption converges to a particular value as κ2 increases,and this asymptote value of κe f f decreases as the switching time increases. These64Figure 5.5: Fitting the probability distribution function in Equation (5.7) toresult calculated from simulation in order to get the effective rate ofabsorption on the boundary.plots suggest that the relationship between the effective rate of absorption, κ2, andswitching time can be approximated by the functionκe f f (α;κ2) = f (α) [1− exp(−g(α)κ2)] , (5.9)where f (α) and g(α) are saturation and rate functions, respectively, to be deter-mined by fitting the model to the curves in Figure 5.6. We call f (α) a saturationfunction because it predicts where the values of κe f f will saturate for each switch-ing time, and g(α) a rate function because it gives the growth rate of κe f f for eachswitching time.Next, we fit the model in Equation (5.9) to each of the curves in Figure 5.6to get the values of the functions f (α) and g(α) for each switching time α . Theresults of the fits are shown in Figure 5.7. We observe from these plots that the65(a) Estimates for half-plane problem. (b) Estimates for disk problem.Figure 5.6: Estimates of the effective rate of absorption κe f f obtained by fit-ting probability distribution functions to the simulated results for differ-ent values of κ2 and switching time.estimates for the functions f (α) and g(α) for the two geometries are very similar.(a) Estimates for half-plane problem. (b) Estimates for disk problem.Figure 5.7: Estimates of the values of the functions f (α) and g(α) obtainedby fitting the model in Equation (5.9) to the curves in Figure 5.6.Lastly, we fit a linear polynomial to the curves of g(α) in Figure 5.7 using the66polyfit command in MATLAB, and for the saturation function, we fit the functionf (α) =µ11+(αµ2) , (5.10)to the curve of f (α) in Figure 5.7 using the least squares curve fitting tool inMATLAB. The parameters µ1 and µ2 are to be determined from the fit.(a) Fits for half-plane problem. (b) Fits for disk problem.Figure 5.8: Fittings functions to the estimated values for the functions f (α)and g(α).For the right half-plane problem, the functions obtained from the fittings aref (α) =0.89701+( α1.7318) , and g(α) = 0.1875α+0.4545 . (5.11)And for the disk problem, we obtainf (α) =0.93441+( α1.6515) , and g(α) = 0.1938α+0.4337 . (5.12)We notice that the functions f (α) and g(α) for the half-plane problem and the diskproblem are very similar. These equations can be used together with the model inEquation (5.9) to predict the effective rate of absorption on the boundary as a func-tion of the rate of absorption κ2 and the switching time, for a situation where theboundary switches between a perfectly reflecting boundary and a specific partially67absorbing boundary.5.2.2 Switching between perfectly reflecting and perfectly absorbingIn this section, we consider a situation where the boundary switches between per-fectly absorbing and perfectly reflecting. Similar to the case of Section 5.2.1, wewant to approximate the switching boundary with a partially absorbing boundary,where the rate of absorption on the boundary is the effective rate of absorption.Following a similar approach to that of Section 5.2.1, we derive a mathematicalrelation for the effective rate of absorption on the boundary as a function of theswitching time. We consider the right half-plane problem and the disk problem,with different switching time. For each switching time, we simulate trajectoriesand calculate the probability of finding a particle performing Brownian motion inthe region of interest. Then we fit the probabilities obtained from the simulation tothe appropriate probability distribution function in order to estimate the effectiverate of absorption on the boundary.(a) Boundary switches every 0.5 time units. (b) Boundary switches every 3 time units.Figure 5.9: Fitting the probability distribution function in Equation (5.7) tothe probability calculated from simulation.Figure 5.9 shows some of the fits of the probability distribution function inEquation (5.7) to simulated results, while Figure 5.10 shows the fits of the proba-bility distribution function in Equation (5.8) to simulated result. For these plots, weaveraged over 10 simulations each with 10,000 particles, D = 0.1, x0 = 5, a = 5,68(a) Boundary switches every 0.5 time units. (b) Boundary switches every 3 time units.Figure 5.10: Fitting the probability distribution function in Equation (5.8) tothe probability calculated from simulation.and r0 = 2.(a) Estimates for half-plane problem. (b) Estimates for disk problem.Figure 5.11: Estimates of the effective rate of absorption κe f f obtained byfitting probability distribution functions to the probability calculatedfrom simulation.For this case, we also considered switching time between 1/2 and 4 with anincrement of 1/2, and for each value of α , we fit the distribution functions to thesimulated probability to get effective rate of absorption κe f f . The results of ourestimates are shown is Figure 5.11. We notice from the curves in this figure that as69the switching time increases, the estimated values of the effective rate of absorptionon the boundary decreases. These plots suggest modelling the relationship betweenthe effective rate of absorption and the switching time with a decaying function.Similar to the results in Section 5.2.1, the estimates of the effective rate ofabsorption for the half-plane problem and the disk problem are very similar. Wealso observe that these curves are similar to the curves for the saturation functionf (α) in Figure 5.7. For this reason, we fit a function similar to the one in Equation(5.10) to these curves;κe f f (α) =γ11+(αγ2) , (5.13)where γ1, and γ2 parameters to be determined from the fit.(a) Fits for half-plane problem. (b) Fits for disk problem.Figure 5.12: Fitting the model in Equation (5.13) to the estimated effectiverates of absorption in Figure 5.12.Fitting the function in Equation (5.13) to the curves in Figure 5.12 using leastsquares curve fitting tool in MATLAB, we obtained the same values for the param-eters γ1 and γ2 for the right half-plane problem and the disk problem. Therefore,the relation for the effective rate of absorption for the two geometry is the sameand it isκe f f (α) =1.00231+( α1.6349) . (5.14)70Thus, we can approximate a boundary that switches between perfectly absorbingand perfecting reflecting with a partially absorbing whose rate of absorption isdetermined using this relation.71Chapter 6Discussion and ConclusionWe have presented a mathematical technique for identifying a partially absorbingbiological boundary, and also techniques for estimating the rate of absorption onsuch a boundary. These techniques are based on the probability of finding a particleperforming Brownian motion in some region of a domain.The first chapter of this thesis gives a brief introduction to biological bound-aries, the kind of biological boundaries we are interested in, and some of the bi-ological terms used in this project. We also presented a brief description of thesingle particle tracking experiment that motivated this project.In the second chapter, we derived the diffusion equation from a simple one di-mensional random walk process and used the derived equation to set-up our mathe-matical problem. Also in this chapter, we gave a brief introduction to a partially ab-sorbing boundary and the Robin boundary condition, and presented some relationsbetween the rate of absorption on the boundary in the Robin boundary conditionand the probability of absorption on the boundary, derived by Erban and Chapman[4], and Steven Andrews [1].The third chapter involves deriving probability distribution functions for find-ing a particle performing Brownian motion in some region of a two dimensionaldomain. We considered an unbounded domain, and a bounded domain in rectangu-lar and polar coordinates. For the unbounded domain, we derived the probabilitydistribution function for finding the particle in some region of the domain. To dothis, we assumed that there is a region of the unbounded domain that is bounded72by an ‘imaginary’ boundary, and then derive the probability distribution functionfor finding the particle in this region at time t. This boundary is imaginary be-cause it does not restrict the motion of the particle in and out of the region. For thebounded domain, we considered the right half-plane with an absorbing boundary atx = 0, and a disk-shaped region also with an absorbing boundary, and then derivedthe probability distribution function for finding the particle in these regions. Sincethese distribution functions give us the likelihood of finding the particle of interestin some region over time, we can fit these functions to experimental data to deter-mine whether or not the experimental particle is experiencing a restricted motion.In order to do this, we need to consider the scenario that best models our exper-imental set-up and collect information such as the initial position of the particle,its diffusion coefficient, and the suspected locations of the boundary from the ex-periment. Substituting these pieces of information into the appropriate distributionfunction gives the distribution of probability for finding the particle in the desig-nated region over time. Using this probability distribution together with the finalposition of the experimental particle, we can conclude whether or not the particleis experiencing a restricted motion and the possibility of a boundary.Supposing we have identified a partially absorbing boundary, the next problemis to estimate the rate of absorption on the boundary. In the fourth chapter, wepresent a technique for estimating the rate of absorption on a partially absorbingboundary. This technique involves deriving a probability distribution function forfinding a particle performing Brownian motion in a region bounded by a partiallyabsorbing boundary. We considered two possible geometries: the right half-planeand a disk-shaped region. The probability distribution functions derived for eachcase contain a parameter that determines the rate of absorption on the boundary.These distribution functions are used together with maximum likelihood estimationto estimate that parameter based on simulated data. To do this, we first simulatedtrajectories using a specified rate of absorption on the boundary, and from thissimulation, we collected the total number of particles that were absorbed on theboundary and those remaining in the region after a specified time. This informationwas substituted into the appropriate likelihood function derived using the relationin Equation (4.39), and then the resulting function was maximized with respect toκ to get the most probable rate of absorption on the boundary. Figures 4.6 - 4.973show the results of our estimates of the rate of absorption on the boundary withstandard deviation error bars for different parameter values. For Figures 4.6 and4.8, we varied the initial position of the particle. We observe from these figuresthat we were able to recover the rate of absorption on the boundary with high levelof accuracy for the cases where the particles start close the boundary, and as theinitial position of the particle moves farther away from the boundary, our estimatesbecome less accurate and the standard deviation increases. This is due to the factthat when the particles start far away from the boundary, they spread out more,therefore, it takes them longer to get to the boundary compared to when they startcloser to the boundary. In Figures 4.7 and 4.9, we fixed the initial position of theparticle and varied the rate of absorption on the boundary. We can see from thesefigures that we have good estimates when the rate of absorption on the boundary issmall, and as κ increase, the accuracy of our estimates decreases and the variancealso increases.We also looked at a scenario where the boundary of the region in which the par-ticle is performing Brownian motion is formed by some other diffusing particles,and as a result of the motion of the particles forming the barrier, the boundary fluc-tuates between a perfectly reflecting boundary and partially absorbing boundaries.To estimate the rate of absorption on the boundary, we used mean field approxima-tion to approximate the fluctuating boundary with a partially absorbing boundarywhere the rate of absorption on the boundary is the average of the fluctuating ratesof absorption. To verify that this approximation works, we simulate trajectoriesin a region bounded by a fluctuating boundary, and calculate the probability thatthe particle performing Brownian motion is still inside the region at time t from thesimulation. In the simulation, the boundary switches from one boundary type to an-other at every unit time interval, and the probability calculated from it is comparedto the result obtained from the probability distribution function derived using themean field approximation (that is, Equation (5.4) for the half-plane problem andEquation (5.6) for the disk problem). The comparison of the two probabilities isshown in Figure 5.2 for the right half-plane problem and in Figure 5.3 for the diskproblem. We notice from these figures that the two result agree to a large extent,this shows that the mean field approximation works for this scenario. Thus, weconclude that the fluctuating boundary in this scenario can be approximated by74a partially absorbing boundary whose rate of absorption on the boundary is theaverage of the fluctuating rates of absorption.Furthermore, we considered a scenario where the boundary switches betweentwo boundary types at specified time interval. The first case considered is when theboundary switches between a perfectly reflecting boundary and a specific partiallyabsorbing boundary (with rate of absorption on the boundary labelled as κ2), whilethe other case is when it switches between a perfectly reflecting and a perfectlyabsorbing boundary. We approximate these switching boundaries with a partiallyabsorbing boundary, where the rate of absorption on the boundary is the rate ofabsorption that best approximates the solution, given the switching on the bound-ary, and it is called the effective rate of absorption denoted by κe f f . For each ofthe two cases, we derived a mathematical relation for finding the effective rate ofabsorption on the boundary a function of the switching time. To get the relationfor the first case, we simulate trajectories in a region where the boundary switchesbetween the two boundary types. From the simulation, we calculate the probabilitythat a particle performing Brownian motion in this region is still in the region attime t, and then fit the appropriate probability distribution function to the calcu-lated probability in order to get the rate of absorption on the boundary that bestfits the simulated result. This rate of absorption is the effective rate of absorptionon the boundary. Figure 5.6 shows the fits for several values of κ2 and differentswitching time, for the half-plane problem and the disk problem. Next, we fit themodel in Equation (5.9) to the estimates in Figure 5.6 in order to get the values forthe saturation function f (α), and the rate function g(α) for each switching time α .And lastly, we fit functions to the curves of f (α) and g(α) obtained from the pre-vious fits. The mathematical relation obtained for the effective rate of absorptionκe f f for the right half-plane problem is given asκe f f (α;κ2) =(0.89701+( α1.7318)) [1− exp(−(0.1875α+0.4545)κ2)] . (6.1)And that of the disk problem isκe f f (α;κ2) =(0.93441+( α1.6515)) [1− exp(−(0.1938α+0.4337)κ2)] , (6.2)75where α is the switching time, and κ2 is the rate of absorption for the partiallyabsorbing boundary that switches.We notice that the relations for the effective rate of absorption on the partiallyabsorbing boundary for the two geometries are very similar. This suggests that therate for this scenario is independent of the geometry of the region, rather it depen-dent on κ2 and switching time. We believe that the slight difference in the results isas a result of approximation errors and the stochastic nature of our approach. Thus,in a situation where a boundary switches between a perfectly reflecting boundaryand a specific partially absorbing boundary, this boundary can be replaced witha partially absorbing boundary where the rate of absorption on the boundary isobtained using the relations in Equations (6.1) and (6.2) appropriately.For the case where the boundary switches between a perfectly reflecting bound-ary and a perfectly absorbing boundary, we used a similar approach of fitting proba-bility distribution functions to the probabilities calculated from simulation in orderto estimate the effective rate of absorption on the boundary. We then fit a functionto these estimates to get a relation that predicts the effective rate of absorption onthe boundary as a function of the switching time. The relation obtained for the righthalf-plane problem and the disk problem are exactly the same for this scenario, andit is given in Equation (5.14). This is in agreement with the suggestion in the pre-vious scenario of switching boundary that the effective rate of absorption on theboundary is independent of the geometry. Therefore, we conclude that a boundarythat switches between a perfectly reflecting and a perfectly absorbing boundaryat regular time intervals can be approximated with a partially absorbing boundarywhose rate of absorption in obtained using the relation in Equation (5.14).We believe that the techniques presented in this project will help in identifyingpartially permeable biological boundaries, and also in estimating the rate of absorp-tion on the boundaries. This will go a long way in enhancing our understandingthe motion of biological particles, and the analysis of experimental data.6.1 Future workFuture work in this direction would be to apply the techniques presented in thisthesis to experimental data.76For the regions we have considered in this thesis, we have assumed that theboundaries are smooth, while in reality, biological boundaries may not be smooth.It would be interesting to develop a technique for estimating the rate of absorptionon a rough boundary.In addition, it would also be interesting to develop a technique for estimat-ing the rate of transmission on a partially transmitting boundary. For a transmittingboundary, a diffusing particle can cross the boundary from both sides of the bound-ary. As a result of this, the boundary has two faces, each with a rate of transmission,and these rates are not necessarily the same. Developing a technique to estimatesthese rates of transmission would also be an interesting question.77Bibliography[1] S. S. Andrews. Accurate particle-based simulation of adsorption, desorptionand partial transmission. Physical Biology, 6(4):046015, 2009. URLhttp://stacks.iop.org/1478-3975/6/i=4/a=046015. → pages 11, 12, 72[2] S. S. Andrews and D. Bray. Stochastic simulation of chemical reactions withspatial resolution and single molecule detail. Physical Biology, 1(3):137,2004. URL http://stacks.iop.org/1478-3975/1/i=3/a=001. → pages 40, 52[3] S. Chandrasekhar. Stochastic problems in physics and astronomy. Reviewsof Modern Physics, 15, 1943. → pages 5[4] R. Erban and S. J. Chapman. Reactive boundary conditions for stochasticsimulations of reaction–diffusion processes. Physical Biology, 4(1):16,2007. URL http://stacks.iop.org/1478-3975/4/i=1/a=003. → pages 11, 32,72[5] F. C. Goodrich. Random walk with semiadsorbing barrier. Journal ofChemical Physics, 22:588–594, 1954. → pages 5[6] G. Lamm. Extended Brownian dynamics. iii. three-dimensional diffusion.The Journal of chemical physics, 80(6):2845–2855, 1984. → pages 5[7] G. Lamm and K. Schulten. Extended Brownian dynamics. ii. reactive,nonlinear diffusion. The Journal of Chemical Physics, 78(5):2713–2734,1983. → pages 5[8] S. D. Lawley, J. C. Mattingly, and M. C. Reed. Stochastic switching ininfinite dimensions with applications to random parabolic PDE. SIAMJournal on Mathematical Analysis, 47(4):3035–3063, 2015. → pages 58[9] N. Ruthardt, D. C. Lamb, and C. Bra¨uchle. Single-particle tracking as aquantitative microscopy-based approach to unravel cell entry mechanisms of78viruses and pharmaceutical nanoparticles. Molecular therapy, 19(7):1199–1211, 2011. → pages 2[10] A. Singer, Z. Schuss, A. Osipov, and D. Holcman. Partially reflecteddiffusion. SIAM Journal on Applied Mathematics, 68(3):844–868, 2008.doi:10.1137/060663258. URL http://dx.doi.org/10.1137/060663258. →pages 5[11] J. St-Pierre and H. L. Ostergaard. A role for the protein tyrosine phosphataseCD45 in macrophage adhesion through the regulation of paxillindegradation. PloS one, 8(7):e71531, 2013. → pages 3[12] Wikipedia. Single particle tracking — wikipedia, the free encyclopedia,2015. URL https://en.wikipedia.org/w/index.php?title=Single particle tracking&oldid=668616350. [Online; accessed 9-May-2016].→ pages 2[13] Wikipedia. PTPRC — wikipedia, the free encyclopedia, 2016. URLhttps://en.wikipedia.org/w/index.php?title=PTPRC&oldid=733512364.[Online; accessed 29-April-2016]. → pages 3[14] Wikipedia. Immunoglobulin g — wikipedia, the free encyclopedia, 2016.URL https://en.wikipedia.org/w/index.php?title=Immunoglobulin G&oldid=731665299.[Online; accessed 29-April-2016]. → pages 3[15] Wikipedia. Likelihood function — wikipedia, the free encyclopedia, 2016.URL https://en.wikipedia.org/w/index.php?title=Likelihood function&oldid=726085070.[Online; accessed 2-May-2016]. → pages 12[16] Wikipedia. Macrophage — wikipedia, the free encyclopedia, 2016. URLhttps://en.wikipedia.org/w/index.php?title=Macrophage&oldid=728586918.[Online; accessed 11-August-2016]. → pages 2[17] Wikipedia. Marcum q-function — wikipedia, the free encyclopedia, 2016.URL https://en.wikipedia.org/w/index.php?title=Marcum Q-function&oldid=717292997.[Online; accessed 3-May-2016]. → pages 19[18] Wikipedia. Random walk — wikipedia, the free encyclopedia, 2016. URLhttps://en.wikipedia.org/w/index.php?title=Random walk&oldid=718681846.[Online; accessed 30-May-2016]. → pages 879[19] E. Zauderer. Partial Differential Equations of Applied Mathematics. JohnWiley & Sons, Inc, CHoboken, New Jersey, Third Edition, 2016. → pages3580Appendix ASupporting MaterialsBelow are some of the smoldyn codes used to produce the results in this thesis.These scripts were called from MATLAB using some other scripts through whichthe value of the parameters such as simulation stopping time, simulation time step,initial number of particles, the staring position of the particles, and the differentrates of absorptions on the boundary were supplied. In addition, the data obtainedfrom these simulations are exported to Matlab for further processing.A.1 Smoldyn code for partially absorbing boundaryproblem# likelihood simulation# rectangular Geometrygraphics opengl#graphic_iter 5dim 2 # specifying the dimension of the systemspecies A # specifying particle species typesdifc all 0.1 # specifying diffusion coefficient of particlecolor A 1 0 0 # specifying the species colordisplay_size all 2.01 # specifying the display size particles81time_start 0 # specifying starting time of simulationtime_stop T_time # specifying stopping time of simulationtime_step T_step # specifying simulation time step### defining system boundaryboundaries 0 0 100 rboundaries 1 0 100 rframe_thickness 0### designing the simulation graphics interfacestart_surface Right_rectangleaction front A reflectcolor back black 0.2color front blue 0.1panel rect +0 0 0 100panel rect -0 100 0 100panel rect +1 0 100 100panel rect -1 0 0 100end_surface### defining the partially absorbing boundarystart_surface Boundaryrate A bsoln fsoln Kappa # specifying the rate of absorption# on the boundaryaction A front reflectcolor back red 0.2panel rect -0 0 0 100thickness 1end_surface### defining simulation region of intereststart_compartment right_compt82surface Right_rectanglepoint 50 50end_compartment### specifying the initial number of particle (called ’total’)### and the start positionmol total A start_point 50### defining the output fileoutput_files FILEROOT_out.txt stdout### writing the total number of molecules left in the region### into the output file at each simulation time stepcmd N 1 molcountincmpts right_compt FILEROOT_out.txtcmd a molcount stdoutend_fileA.2 Smoldyn code for fluctuating boundary problem# Fluctuating boundary simulation# rectangular Geometrygraphics opengl#graphic_iter 5dim 2 # specifying the dimension of the systemspecies A # specifying particle species typesdifc all 0.1 # specifying diffusion coefficient of particlecolor A 1 0 0 # specifying the species colordisplay_size all 2.01 # specifying the display size particlestime_start 0 # specifying starting time of simulation83time_stop T_time # specifying stopping time of simulationtime_step T_step # specifying simulation time step### defining system boundaryboundaries 0 0 100 rboundaries 1 0 100 rframe_thickness 0### designing the simulation graphics interfacestart_surface Right_rectangleaction front A reflectcolor back black 0.2color front blue 0.1panel rect +0 0 0 100 r1panel rect -0 100 0 100 r2panel rect +1 0 100 100 r3panel rect -1 0 0 100 r4end_surface### defining the boundary that fluctuatesstart_surface Boundaryaction A front reflectcolor back red 0.2panel rect -0 0 0 100 rec1thickness 1end_surface## specifying an initial rate of absoprtion on the boundarysurface Boundary rate A bsoln fsoln 0### defining simulation region of intereststart_compartment right_comptsurface Right_rectangle84point 50 50end_compartment### specifying the initial number of particle (called ’total’)### and the start positionmol total A start_point 50### defining the output fileoutput_files FILEROOT_out.txt stdout### performing the switching on the boundarycmd @ 1 set surface Boundary rate A bsoln fsoln aacmd @ 2 set surface Boundary rate A bsoln fsoln bbcmd @ 3 set surface Boundary rate A bsoln fsoln cccmd @ 4 set surface Boundary rate A bsoln fsoln ddcmd @ 5 set surface Boundary rate A bsoln fsoln ee: : : : : : : : : : : : : : : : : : : : :. . . . . . . . . . . . . . . . . . . . .cmd @ 76 set surface Boundary rate A bsoln fsoln a75cmd @ 77 set surface Boundary rate A bsoln fsoln a76cmd @ 78 set surface Boundary rate A bsoln fsoln a77cmd @ 79 set surface Boundary rate A bsoln fsoln a78cmd @ 80 set surface Boundary rate A bsoln fsoln a79### writing the total number of molecules left in the region### into the output file at each simulation time stepcmd N 1 molcountincmpts right_compt FILEROOT_out.txtend_file85
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Mathematical modelling of partially absorbing boundaries...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Mathematical modelling of partially absorbing boundaries in biological systems Iyaniwura, Sarafa Adewale 2016
pdf
Page Metadata
Item Metadata
Title | Mathematical modelling of partially absorbing boundaries in biological systems |
Creator |
Iyaniwura, Sarafa Adewale |
Publisher | University of British Columbia |
Date Issued | 2016 |
Description | This project presents a mathematical framework for identifying partially permeable biological boundaries, and estimating the rate of absorption of diffusing objects at such a boundary based on limited experimental data. We used partial differential equations (PDEs) to derive probability distribution functions for finding a particle performing Brownian motion in a region. These distribution functions can be fit to data to infer the existence of a boundary. We also used the probability distribution functions together with maximum likelihood estimation to estimate the rate of absorption of objects at the boundaries, based on simulated data. Furthermore, we consider a switching boundary and provide a technique for approximating the boundary with a partially permeable boundary. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2016-08-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0308701 |
URI | http://hdl.handle.net/2429/58907 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of |
Degree Grantor | University of British Columbia |
GraduationDate | 2016-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2016_september_iyaniwura_sarafa.pdf [ 1.81MB ]
- Metadata
- JSON: 24-1.0308701.json
- JSON-LD: 24-1.0308701-ld.json
- RDF/XML (Pretty): 24-1.0308701-rdf.xml
- RDF/JSON: 24-1.0308701-rdf.json
- Turtle: 24-1.0308701-turtle.txt
- N-Triples: 24-1.0308701-rdf-ntriples.txt
- Original Record: 24-1.0308701-source.json
- Full Text
- 24-1.0308701-fulltext.txt
- Citation
- 24-1.0308701.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0308701/manifest