Mathematical models of motor-based intracellular transport

The efficient transport of particles throughout a cell plays a fundamental role in several cellular processes. Broadly speaking, intracellular transport can be divided into two categories: passive and active transport. Whereas passive transport generally occurs via diffusive processes, active transport requires cellular energy through adenosine triphosphate (ATP). Many active transport processes are driven by molecular motors such as kinesin and dynein, which carry cargo and travel along the microtubules of a cell to deliver specific material to specific locations. Breakdown of molecular motor delivery is correlated with the onset of several diseases, such as Alzheimer's and Parkinson's. We mathematically model two fundamental cellular processes. In the first part, we introduce a possible biophysical mechanism by which cells attain uniformity in vesicle density throughout their body. We do this by modeling bulk motor density dynamics using partial differential equations derived from microscopic descriptions of individual motor-cargo complex dynamics. We then consider the cases where delivery of cargo to cellular targets is (i) irreversible and (ii) reversible. This problem is studied on the semi-infinite interval, disk, and spherical domains. We also consider the case where exclusion effects come into play. In all cases, we find that allowing for reversibility in cargo delivery to cellular targets allows for more uniform vesicle distribution. In the second part, we see how active transport by molecular motors allows for length control and sensing in flagella and axons, respectively. For the flagellum, we model length control using a doubly stochastic Poisson model. For axons, we model bulk motor dynamics by partial differential equations, and show how spatial information may be encoded in the frequency of an oscillating chemical signal being carried by dynein motors. Furthermore, we discuss how frequency-encoded signals may be decoded by cells, and how these mechanisms break down in the face of noise.

MATHEMATICAL MODELS OF MOTOR-BASED INTRACELLULAR TRANSPORT by Bhargav Karamched A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Mathematics The University of Utah May 2017 c Bhargav Karamched 2017 Copyright All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of Bhargav Karamched has been approved by the following supervisory committee members: Paul Bressloff , Chair(s) 27 February 2017 Date Approved James Keener , Member 27 February 2017 Date Approved Aaron Fogelson , Member 27 February 2017 Date Approved Fernando Guevara-Vasquez , Member 27 February 2017 Date Approved Michael Vershinin , Member 27 February 2017 Date Approved by Peter Trapa , Chair/Dean of the Department/College/School of Mathematics and by David B. Kieda , Dean of The Graduate School. ABSTRACT The efficient transport of particles throughout a cell plays a fundamental role in several cellular processes. Broadly speaking, intracellular transport can be divided into two categories: passive and active transport. Whereas passive transport generally occurs via diffusive processes, active transport requires cellular energy through adenosine triphosphate (ATP). Many active transport processes are driven by molecular motors such as kinesin and dynein, which carry cargo and travel along the microtubules of a cell to deliver specific material to specific locations. Breakdown of molecular motor delivery is correlated with the onset of several diseases, such as Alzheimer's and Parkinson's. We mathematically model two fundamental cellular processes. In the first part, we introduce a possible biophysical mechanism by which cells attain uniformity in vesicle density throughout their body. We do this by modeling bulk motor density dynamics using partial differential equations derived from microscopic descriptions of individual motor-cargo complex dynamics. We then consider the cases where delivery of cargo to cellular targets is (i) irreversible and (ii) reversible. This problem is studied on the semiinfinite interval, disk, and spherical domains. We also consider the case where exclusion effects come into play. In all cases, we find that allowing for reversibility in cargo delivery to cellular targets allows for more uniform vesicle distribution. In the second part, we see how active transport by molecular motors allows for length control and sensing in flagella and axons, respectively. For the flagellum, we model length control using a doubly stochastic Poisson model. For axons, we model bulk motor dynamics by partial differential equations, and show how spatial information may be encoded in the frequency of an oscillating chemical signal being carried by dynein motors. Furthermore, we discuss how frequency-encoded signals may be decoded by cells, and how these mechanisms break down in the face of noise. CONTENTS ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii NOTATION AND SYMBOLS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii CHAPTERS 1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Microtubules, Molecular Motors, and Intracellular Transport . . . . . . . . . . . . . 2 1.1.1 Cytoskeleton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Kinesin and Dynein . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.3 Neurons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Piecewise Deterministic Markov Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2.1 Adiabatic Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.3 Chemical Master Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.4 Structure of Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2. CELL GEOMETRY AND VESICULAR DELIVERY . . . . . . . . . . . . . . . . . . . . . . . . 16 2.1 Semi-infinite Track . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Irreversible Delivery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Reversible Delivery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Cayley Tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Irreversible Delivery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Reversible Delivery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Higher-dimensional Geometries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 The Disk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 The Sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Discrete Microtubule Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 The Disk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 The Sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. 18 18 19 21 23 25 27 28 31 34 34 38 42 VESICLE DELIVERY WITH EXCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.1 Irreversible Vesicular Transport with Exclusion . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Steady-state Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Reversible Vesicular Transport with Exclusion . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Mean-field and Continuum Limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Fast Switching Limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Steady-state Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.4 Disparity in Hopping Rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 49 51 54 55 58 62 3.3 Relationship with Other Exclusion Process Models . . . . . . . . . . . . . . . . . . . . . . 63 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4. FLAGELLAR LENGTH CONTROL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.1 Stochastic Model of Flagellar Length Control . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 IFT Injection as a Poisson Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 Stochastic Model of IFT Docking at the Basal Body . . . . . . . . . . . . . . . . . 4.1.3 RanGTP Model of IFT Flux Regulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Analysis of Model Using the Theory of Doubly Stochastic Poisson Processes 4.2.1 Doubly Stochastic Poisson Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Calculation of Mean and Variance of IFT Numbers within Flagellum . . 4.2.3 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Modified Stochastic Model for Small M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Feynman-Kac Formula and Reduction to a DSPP . . . . . . . . . . . . . . . . . . 4.3.2 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. 72 74 75 77 79 79 82 84 86 87 89 90 AXONAL LENGTH SENSING . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.1 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.1.1 Delayed Feedback Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.1.2 Linear Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.1.3 Advection-Diffusion Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2 Analysis of PDE Model Using Green's Functions . . . . . . . . . . . . . . . . . . . . . . . 102 5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.3.1 Knockdown of Molecular Motor Activity . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.3.2 Effects of Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.3.3 Numerical Methods and Parameter Values . . . . . . . . . . . . . . . . . . . . . . . . 110 5.4 Frequency Decoding by a Feedforward Gene Network . . . . . . . . . . . . . . . . . . 112 5.5 Effects of Intrinsic Noise on Axonal Length Sensing . . . . . . . . . . . . . . . . . . . . . 116 5.5.1 Errors in Axonal Length Sensing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6. FUTURE DIRECTIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 6.1 6.2 6.3 6.4 Microscopic Descriptions of Motor Transport . . . . . . . . . . . . . . . . . . . . . . . . . . 124 Motor Routing by Microtubule Network Modification . . . . . . . . . . . . . . . . . . . 125 Neuron Polarization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 Closing Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 APPENDIX: GREEN'S FUNCTION FOR ADVECTION DIFFUSION EQUATION . 129 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 SUBJECT INDEX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 v NOTATION AND SYMBOLS R Set of real numbers Rn n-dimensional Cartesian space df dy Derivative of the function f with respect to the variable y ∂f ∂y Partial derivative of the function f with respect to the variable y |x| absolute value of a number x ∈ R ||x|| Euclidean norm of a vector x ∈ Rn h f , gi inner product of functions f , g, h f , gi ≡ fˆ Fourier transform of a function f f˜ Laplace transform of a function f ( f ∗ g)( x ) Convolution of two functions, ( f ∗ g)( x ) ≡ Br (σ) ball in Rn centered at σ with radius r > 0 h N i, E[ N ] expected value of some random variable N h MN i covariance of two random variables M, N Var[ N ] variance of the random variable N σN standard deviation of random variable N, σN ≡ CV coefficient of variation of a random variable, CV ≡ exp( x ) natural exponential function, exp( x ) ≡ e x Kν ( x ) Modified Bessel function of the second kind of order ν δ( x ) Dirac-Delta function H (x) Heaviside function R f ( x ) g( x )dx R f (y) g( x − y)dy p Var[ N ] σN hNi ACKNOWLEDGEMENTS I've been blessed to be part of a math department whose members are bright and supportive and to be from a family that pushed me to be my best. This dissertation may have my name on it, but in actuality, it is the fruit of a sequence of collaborations with several different people. Here is an effort to express my gratitude to several people (though not all!) who have helped me reach this point in my life. First of all, thanks to my parents, whose fortitude in the face of several crises inspired me on several occasions to get back up after each and every failure. Amma, thank you for teaching me that failure is nothing to be afraid of, so long as you try your best. Nanna, thank you for making me laugh and see the goodness in life even during my darkest times. Thanks to my sister, Keerthi, for your love and your eternal support for my work. Thank you for being excited when I told you about my results even though you don't understand anything about my results. None of this would have been possible without my friends' help. Thank you Leif, for being a great office mate, distracting me from the difficulties of my research, and engaging me in conversations about sports and obscure fields of mathematics. These three years have been the most productive of my life, and a lot of it is due to you. Thank you Shiu-Tang for giving life to the adventurer inside of me and taking me on hikes in some of the most beautiful places on Earth. Thank you Pavel for showing me how to laugh at the absurdity of life. Our workouts together have taught me that I'm stronger than I give myself credit for, physically and emotionally. Utah has a wonderful faculty. Thanks to Prof. Jim Keener for recruiting me, introducing me to the wonderful field of Mathematical Biology, and giving me insights into mathematical problems when I walked into your office with them. You've shown me that true passion for a subject need not diminish as time progresses. Thanks to Prof. Fernando Guevara-Vasquez for giving me great encouragement and advice during my first year when I TA'd for you. Thanks to Will Nesse for giving me the opportunity to teach great classes. Last, but most definitely not least, thanks to my advisor Prof. Paul Bressloff. You've taken this shy, unsure Indian boy and guided him to a place he's only dreamed of since he was 10 years old. Thanks for your patience with my constant visits to your office and my inconsistency with parameter values. You're such an inspiration. Seeing you and your work made me happy about my decision to come to Utah. You've made me see the potential I have and what I am capable of doing with it. I can only hope to do for another student in the future what you have done for me. Thanks for the coffee. You're the best. viii CHAPTER 1 INTRODUCTION One of the major challenges of modern biology is to understand molecular motor-based intracellular transport and the biophysical processes that are facilitated by it. This involves many different levels of description, from the biochemical reactions that dictate an individual motor's movement and vesicular uptake to bulk molecular motor movement across entire cells for delivery of these vesicles to specific locations within the cell. From a mathematical perspective, motor-based intracellular transport provides a rich area for application of ordinary, delay, and partial differential equations, statistical physics, and the theory of stochastic processes. The latter are particularly pertinent due to the interior of a cell being a crowded, fluctuating, anisotropic environment [18, 95]. Interest in intracellular transport processes has steadily grown in a diversity of academic fields over the years (see, for example, [3, 11, 14-18, 21, 46, 56, 57, 59, 63, 65, 72, 76, 78, 81, 83, 87, 89, 102, 125, 127]). In this dissertation, we introduce mathematical models that address two fundamental biological questions that have been experimentally observed to be intertwined with motorbased intracellular transport processes [3, 78, 102, 125]: 1. How are cells able to route motors so that they deliver specific cargo to particular subcellular compartments in the face of molecular crowding and anisotropy? In cases where vesicular density throughout a cell needs to be relatively uniform, how is this homogeneity achieved? 2. How do cells that are relatively large (e.g., neurons) control and sense their own size in the face of noisy intracellular environments? We note that this work focuses on bulk motor travel along microtubules and not on the biochemical reactions that modulate an individual motor's movement and uptake of vesicles. 2 In order to put our mathematical models into context, we will provide a brief background of the biological processes and structures that form the basis for the phenomena we model in the following. Later, we will describe the fundamental mathematical structures underlying several of the equations we use for modeling bulk motor motion. 1.1 Microtubules, Molecular Motors, and Intracellular Transport Broadly speaking, there are two basic mechanisms for intracellular transport: passive diffusion within the cytosol or the surrounding plasma membrane of the cell and active motor-driven transport along polymerized filaments such as microtubules and F-actin that comprise the cytoskeleton [18]. Active transport is necessary for the movement of material across relatively long distances, as a diffusion-based mechanism for such transport would require a significantly longer time to reach a particular distance than a given cell allots [60]. This problem becomes particularly acute for large, highly branched cells such as neurons, which range from a micron to a meter in length in humans. Hence, active transport is essential for the efficient delivery of proteins and molecular products to their correct location within a given cell, which is necessary for healthy cellular function and development [1]. Indeed, breakdown in intracellular active transport has been implicated in diseases such as Alzheimer's, Parkinson's, and others [123]. There are two central molecular players in active intracellular transport: (i) the cytoskeleton and (ii) molecular motors. 1.1.1 Cytoskeleton Most proteins are globular proteins, meaning they harvest and store free energy, transform biological compounds into others, or decode genetic information. But many of the most abundant proteins are fibrous proteins that are elongated and often insoluble; these proteins determine the shape and other physical attributes of cells and organisms. Collectively, these proteins that form a cell's structural components are called the cytoskeleton. A typical eukaryotic cell contains three types of cytoskeletal proteins that form fibers extending throughout the cell: microfilaments (with diameter of about 0.7 nm), intermediate filaments (10 nm), and microtubules (24 nm). Although all these proteins provide structural integrity for a cell, their microscopic configurations are markedly different, causing them 3 to have distinct physical properties. Microfilaments, for example, are composed of actin monomers whereas microtubules are composed of tubulin dimers. The term "intermediate refers to an assortment of protein fibers that form subset 16filaments" 8 4 Polaym ers andofMthe oleccytoskeleton ular Motors whose range lies between 0.7 and 24 nm. They are composed of a variety of proteins [95]. left-half complex plane. Finally, DJ has a center at λ = − 1 and radius r, which also One important characteristic of microfilaments and microtubules that differentiates them lies the left-half complex plane provided that r ≤ 1. ± from intermediate filaments is their participation in unbind the dynamic functions of a rates givenkcell Now suppose that actin monomers can bind or at both ends with on ± and , asfilaments shown inexhibit Fig. 4.5. The meaning binding rate is multiplied by acan fixed [43].koff Both polarity, one end of the filament be background distinguished monomer concentration a. (The spatial effects of a nonuniform monomer concenfrom the other based on the orientation of the individual components building up the tration are considered by Edelstein-Keshet and Ermentrout [157]; see also Ex. 4.3.) filament. Actin monomers have oneends end with an to ATP-binding cleft. The microfilament end The difference between the two is due the fact the ATP-actin quickly hydrolyzes ATD-actin so that the consists ofother ATP-actin andThe thebuilding tail consists exposingto these clefts is designated as tip the (-) end and end as (+). blocks of ATD-actin. Rather than writing down the master equation for the system, let us of microtubules are tubulin dimers, one monomer beingnα-tubulin and the other consider the equations for the meanwith number of monomers ± added at each end. Assuming is sufficiently long, we have is designated as the (+) end, β-tubulin. that Thethe endfilament of the microtubule exposing α-tubulin and the other is the (-)dn end. The (+) and (-) designations in each case refer to the tendency dn + − = k+ a − k+ , = k− a − k− . (4.1.9) on of the filament to synthesize ononthe (+)offend and at off the (-) end, although both dt dt dissociate and dissociation occur on both ends; see Figure 1.1 and Figure −1.2. Microtubules − . If a+ ≈ a±c = koff / kon Itsynthesis is clear that the ± end grows provided that a > a±c , where c − , then + − microfilaments are dynamic structures, meaning they are constantly growing aand both ends shrink or grow simultaneously. On the other hand, if a < a < aand c c c then the plus end grows at the same time the minus end shrinks. Finally, adding the dissociating, and when the rate of synthesis at the (+) end of a polarized filament is equal pair of Eq. (4.1.9) shows that to the rate of dissociation at the (-) end, the filament is said to be treadmilling. Hence, a given cell can quickly alter rates ofdn synthesis in response to some external = k aor − dissociation k , on dt off stimulus to modify the length of a filament. In this manner, a cell may use these filaments + − + − . Hence, if the monomer n− , koffMicrofilaments = koff + koff , are andbetter kon =forkon kon with n = n+for+motility. as a means cell+motility, as they are like thin rods concentration a = a0 , where and malleable. Microtubules are hollow tube-like structures and therefore more resistant + − koff + koff structural koff to bending; they are better suited for maintaining integrity; see Figure 1.2. Their a0 = , + − k k + k on on on high abundance, large surface area, and rigidity make them an ideal means of transport- _ + Fig. 4.5: Model of F-actin undergoing polymerization at both ends Figure 1.1. Figure depicting polymerization and dissociation at both ends of a polarized filament. Adapted from [12]. then the total filament length remains constant even though monomers are constantly moving along its length-treadmilling. found as large 13-stranded (each strand is called a protofilament) hollow tube structures. Also, the a and b tubulin used for building the microtubules not only alternate, but they are actually added in pairs. Both the a-tubulin and b-tubulin must bind to GTP to associate, but once bound, the GTP bound to a-tubulin does not move. On the other hand, GTP bound in the b-tubulin may be hydrolyzed to GDP. GDP-bound ab-dimers will not be added to a microtubule, so similar to the situation with ATP and g-actin, if the tubulin has GDP bound to it, it must first exchange it for a GTP before it can be polymerized. Although the affinity of tubulin for GTP is higher than the affinity for GDP, this process is usually facilitated by a GEF, or guanine nucleotide exchange factor. As the signal transduction chapter will show in more detail, this type of nucleotide exchange is a common mechanism for activation of various biochemical pathways. = Microtubules Figure 4. Microtubules exhibit dynamic instability. GTP-bound ab-tubulin dimers are added onto the microtubule. Once the GTP is hydrolyzed, the conformational shift strains the microtubule, which will tend to break apart unless new tubulin dimers are added to stabilize the structure. = Intermediate Filaments = Actin Filaments GTP 4 GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP GTP P GDPGT GTP GTP GDP GDP GDP GDPGDP GDP GDP GDP GDP GDP GDP GDPGDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDPGDP GDP GDP GDP GDP GDP GDP GDPGDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GTP GTP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP P GDP GD GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP P GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GD GD GDP GDP GD GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP P GDP P GDP GDP GDP GDP GDP GDP GD P GDP GD P GDP GDP GDP GDP P GD GDP GDP GDPGDPGDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GDP GD P like nucleus. actin, the tubulin keletal element distribution in a prototypical eukaryotic cell. The purple ballAgain is the itself has enzymatic activity, and over time, the GTPase activity hydrolyzes the GTP to GDP and phosphate. This changes the attachment between b-tubulin of one dimer and the a-tubulin of the dimer it is stacked on because the shape of the subunit changes. Even though it isn't directly loosening its hold on the neighboring tubulin, the shape change causes increased stress as that part of the microtubule tries to push outward. This is the basis of a property of microtubules known as dynamic instability. If there is nothing to stabilize the microtubule, large portions of it will fall apart. However, as long as new tubulin (which will have GTP bound) is being added at a high enough rate to keep a section of low-stress "stable"-conformation microtubule (called the GTP cap) on top of the older GDP-containing part, then it stabilizes the overall microtubule. When new tubulin addition slows down, and there is only a very small or nonexistent cap, then the microtubule undergoes a catastrophe Figure 1.2. Figure depicting structure of microtubules and example of orientation in cells. ajor components of the cytoskeleton are microtubules, microfilaments, and Public downloaded from Wikipedia Commons. e filaments. Each of these are domain polymersfigure composed of repeating subunits in ngements. With just a quick glance (fig. 1), it is very clear that the interments will likely play a significantly different role from either microtubules ing throughout a cell. that They act as a network of roads for the cell ments. Because the IF's arematerial made of long fibrous subunits coileffectively around to form the filament,[95]. there Indeed, is clearly microtubules a great deal of contact interact(which with facilimolecular motors such as kinesin and dynein ion of hydrogen bonds, aka molecular velcro™) between subunits providChapter 12, Cytoskeleton, version 1.0 todifficult facilitate transport of vesicles across relatively large distances. Hence, microtubules nsile strength. It is very to break these subunits apart, and thus the marily used for long-term or permanent load-bearing purposes. Looking at an integral in active intracellular Throughout this dissertation, when wo components of theplay cytoskeleton, one role can see that with the globulartransport. inous shape of the subunits, the maximum area ofdomain contact between subunits we describe a spatial upon which bulk motor dynamics are modeled, we are mited (think of the contact area when you push two basketballs together), asier to separate the subunits or break the microfilament microtubule. effectively abstracting the space or covered by microtubules and giving it a mathematical n use this characteristic to its advantage, by utilizing these kinds of cyformalism. bers in dynamic situations where formation or destruction of intermediate ould take far too long. We now address these three groups of cytoskeletal Microtubules in a given animal cell originate at an organelle near the cell's nucleus more detail. te Filaments known as a centrosome, which contains the microtubule organizing center (MTOC). They emanate from this organelle outward toward the cell membrane. Generally, they are ori- te filaments" is actually a generic name for(+) a family of proteins ented so that their ends are nearer (grouped the cell membrane. In highly polarized cells such es based on sequence and biochemical structure) that serve similar funcMost intermediate filaments fall between 50-100 kDa, including ascell neurons, microtubules are generally oriented that their ends are (60-70 at thekDa), axonal ecting and shaping the or its components. Interestingly, they can even so keratins (40-67(+) kDa), lamins and neurofilaments (62ide the nucleus. The nuclear lamins, which constitute class V intermediate 110 kDa). Nestin (class VI), found mostly in neurons, is an excepterminus. The collective set of microtubules form ation, dense, complex network, and may be orm a strong protective mesh attached to the inside face of the nuclear at approximately 240 kDa. ytoskeleton, version 1.0 navigated by molecular motors such as kinesin and dynein to deliver vesicles to specific Page 176 intracellular locales; see Figure 1.3. 1.1.2 Kinesin and Dynein Kinesin and dynein, proteins that navigate microtubule networks to transport vesicles across cells, are examples of a wider class of proteins called molecular motors. They are biological machines that facilitate any type of movement in living organisms [1, 95]. Examples include myosins, which aide in muscle contraction, topoisomerase, which reduces 5 is, and are Figure 1.3. Micrograph of nucleus, MTOC, and microtubule network. Public domain figure downloaded from Wikipedia Commons. This photo shows the microtubules in a cell. The MTOC (microtubule organising centre) can be supercoiling of DNA in cells, DNA/RNA polymerase, which reads DNA and creates a second strand or transcribes RNA, and kinesin and dynein. There is a vast and rich litera- seen. bar, interested 20 µm).in Another ture on these proteins. We(Scale are specifically incorporatingfluorescent kinesin and dynein proteins in our mathematical models because of their integral role in active intracellular transport. picture of microtubules, which contains dividing Kinesin is a relatively cells.large protein, with a molecular mass of 380 kD. It has two large globular heads and coiled tail domain; see Figure 1.4. Each 10-nm-long head contains a tubulin-binding site and a nucleotide-binding site. The tips of the coiled tails bind to proteins in the membranes of a vesicle, which is a small structure inside of the cell ade up of tubulin dimers as shown in the diagram (left). The molecule of alpha-tubulin and one of beta-tubulin. Each that bindbind the motor on one and aand cargo molecule or or les that the motor on side, one side, a cargo molecule ther examination the cargo androuting the routing the by cargo by examination of theofcargo and the of the of cargo )s was discussed the vesicular transport chapter. discussed in theinvesicular transport chapter. 6 n (B) are are g microcro move to to ethe to the requires res sites are are AA Heavy Chain ATP Heavy Chain ATP Light Chain Light Chain ATP ATP ATP B B ATP Heavy Chain Heavy Chain ATP ATP Intermediate Chain Intermediate Chain ow) by yue) ide the ATP and molin Light Chain Light Chain Figure 1.4. Kinesin (top) and Dynein (bottom). Redrawn from [95]. consisting of fluid enclosed in a lipid bilayer. They can be thought of as packages filled with contents to be shipped to a given locale. The vesicle and its contents become kinesin's cargo [95]. There are several proposed mechanisms for how kinesin moves along microtubules + - [120, 124], but there is one that is widely accepted, which we describe in the following. _ _ + Kinesin moves along a microtubule utilizing a mechanism that involves coordination between its two large globular heads. One head will bind to a β-tubulin with its tubulinbinding site and then adenosine triphosphate (ATP) will bind to the tubulin-bound head. This biochemical binding induces a conformational change in the kinesin and thrusts the unbound head forward; the process then repeats. Hence, kinesin literally walks along microtubules to transport cargo [1, 95]. Dynein has a similar structure to kinesin, but is significantly larger, with a molecular mass of approximately 1.5 MD. It is a dimer of dimers. Two identical heavy chains contain a nucleotide-binding sites, and are responsible for harnessing free energy to change protein Type I myosin, priA B conformation and propel the motor forward. Light chains on the opposite end of the o f-actin, including protein bind to the proteins in the membrane of a vesicle, and are responsible for keeping yosin, binds I myosin, prif-actin A B against each other. tin, including C ular transport. binds f-actin (D) aytosis. i n st ea(E) ch oType t her .XI ytoplasmic stream- 7 dynein attached to its cargo [64, 121]; see Figure 1.4. It is largely unclear what the mechanism for the cooperation between the heavy chains is to propel dynein forward, but it is assumed to be similar to kinesin's mechanism described above, with the underlying factor being the bending energy associated with ATP-induced conformational change [119]. Kinesin tends to walk toward the (+) end of a microtubule whereas dynein walks toward the (-) end [48]. In the context of generic cells, this means kinesin tends to carry cargo toward the cell membrane and dynein tends to bring cargo toward the nucleus. In axonal transport, kinesin tends to haul cargo in the anterograde direction; dynein moves cargo in the retrograde direction. Both motors dynamically bind and unbind completely from a given microtubule; it is relatively uncommon for a motor protein to stay put on a single microtubule for an extended period of time. Furthermore, several dynein and kinesin motors can bind to a single cargo element at the same time, and as these motors pull in opposite directions, the result is a tug-of-war battle to move a cargo in a given direction [63, 87]. It is therefore typical to observe a cargo element undergoing biased bidirectional motion on a given microtubule. Indeed, experimental observations of flourescently tagged cargo elements in nerve cells of Drosophila and C. elegans reveal that cargo elements exhibit ballistic motion in either the retrograde or anterograde directions interspersed with long time periods where cargo elements are stationary; see Figure 1.5 [78, 125]. These experiments show that anterograde and retrograde cargo velocities are on the order of 1 µm/sec. Hence, velocities in all our models of motor transport are taken to be 1 µm/sec. 1.1.3 Neurons Although the biology described above is true for general cells, we have taken care to bring up neurons as a means of context. This is because several of the models we have developed describe molecular motor dynamics and transport processes in neurons. All problems associated with transport processes in cells are even more acute for neurons due to their size and complex, branched structure. Hence, mathematical models of transport processes in neurons lead to very interesting mathematics; they are especially emphasized throughout this thesis, although others are addressed in Chapter 2 and Chapter 4. 8 Differential Regu Figure 1.5. Kymograph of cargo in axon of C. elegans, showing cargo undergo bidirectional motion. Bar graph shows velocity of cargo is on the order of 1 µm/sec. Adapted from [78]. 9 1.2 Piecewise Deterministic Markov Process When considering the active transport of intracellular cargo over relatively long distances (i.e., axons), it is convenient to ignore the microscopic details of how a motor performs a single step and to instead focus on the transitions between states of anterograde and retrograde transport or outward and inward transport. This has motivated a class of macroscopic models which capture the essence of cargo motility as described in [78], wherein there are 3 states of a particle: (1) particle moves ballistically in anterograde direction, (2) particle moves ballistically in retrograde direction, and (3) particle is stationary; see Figure 1.6. In each state, dynamics of the particle are deterministic, but transitions between each state are governed by a Markov process. A system exhibiting the aforementioned properties is called a piecewise deterministic Markov process or a stochastic hybrid system. Consider a particle moving according to the three-state model described above on a one-dimensional (1D) track of length L. As mentioned before, this domain is interpreted as representing a microtubule. Within the track, 0 < x < L, the particle is taken to be in one of three states labeled by n = 0, ±: stationary, moving to the right (anterograde; outward) with speed v+ , or moving to the left (retrograde; inward) with speed v− † . For simplicity, assume v+ = v− , which is justified by the findings in [78]. Let X (t) and N (t) denote the random position and state of the particle at time t and define P( x, n, t|y, m, 0)dx as the joint probability that x < X (t) < x + dx and N (t) = n given that initially the particle was at random position X (0) = y and in the state N (0) = m. Setting pn ( x, t) ≡ ∑ P(x, t, n|0, 0, m)σm , (1.1) m with initial condition pn ( x, 0) = δ( x )σn , ∑m σm = 1, the evolution of the probability is described by the following piecewise deterministic Markov process (PDMP) for t > 0: ∂p+ ∂p+ = −v − β + p+ + αp0 , ∂t ∂x ∂p− ∂p− =v − β − p− + αp0 , ∂t ∂x ∂p0 = β + p+ + β − p− − 2αp0 . ∂t (1.2) (1.3) (1.4) † Due to our emphasis on axonal transport, we will generally call the (+) direction anterograde and the (-) direction retrograde. mobile states, andD 0 is the diffusivity in the unbound state Substituting such a solution into Eqs.(1) and performing n ¼ 0. We are also assuming that there is a uniform, an asymptotic expansion inwn then leads to Eq.(2) to continuous distribution of presynaptic targets along a leading order in . In particular, D − D 0 ρ0 ¼ Oð Þ. region A of the axon, and that the motor complex can Now, suppose that we have a population of motor irreversibly deliver its SVP to a presynaptic target at a complexes injected at one end of the 10 axon at a fixed rate uniform rateκ. Thus, χ A ðxÞ denotes an index function with J 1 , each of which carries a single SVP. As a further microtubule _ + motor complex n=0 α β α β n=− v _ n=+ v _ + + FIG. 1 (color online). Three-state model of the bidirectional transport of a single motor-cargo complex. The particle switches between ðn ¼ −Þof speedv − . The an anterograde state ðn ¼ þÞ of speedv þ , a stationary or slowly diffusing staten( ¼ 0), and a retrograde state motor complex can only deliver a SVP to a presynaptic target in the state n ¼ 0. Figure 1.6. Illustration of three-state model. Adapted from [12]. 168101-2 The parameters α, β ± are the transition rates between the stationary and motile states. We supplement the above PDMP with appropriate boundary conditions at x = 0, L. For example, a reflecting condition at x = 0 and absorbing condition at x = L means that p− (0, t) = p+ (0, t), p− ( L, t) = 0. (1.5) In the general case where v+ 6= v− , transport will be biased toward the right (left) if v+ /β + > v− /β − (v+ /β + < v− /β − ). We note that the above PDMP can be written compactly in matrix-vector format: ∂p = L(p) + Ap, ∂t with   p+ p ≡  p−  , p0   −β+ 0 α −β− α , A≡ 0 β+ β − −2α (1.6)   −v ∂∂xf1   L(f) ≡  v ∂∂xf2  . 0 (1.7) The three-state model is the microscopic description of individual motor-cargo dynamics that underlies several of the partial differential equation (PDE) models we use in our mathematical models. However, we take a moment here to point out that the three-state model is a special case of a general class of motor transport models, in which there are N distinct velocity states labeled n = 1, . . . , N with corresponding velocities vn . For the sake of illustration, we will show how a PDMP model may be used to describe tug-of-war dynamics described in the previous section. 11 Consider a motor complex consisting of N+ anterograde motors and N− retrograde motors. At a given time t, the internal state of the motor-cargo complex is completely characterized by the ordered pair (n+ , n− ) of the numbers of anterograde and retrograde motors that are bound to the microtubule and therefore actively pulling the cargo. The velocity v(n+ , n− ) of a given state is given by [12, 18] v(n+ , n− ) = n+ Fs+ − n− Fs− , n+ Fs+ /v f + + n− Fs− /vb− (1.8) where Fs± is the stall force for anterograde and retrograde motors, i.e., the magnitude of force opposing a single motor's movement that renders the single motor's velocity 0, v f + is the anterograde motor's velocity in the absence of any opposing force and vb− is the retrograde motor's velocity in the anterograde direction when the opposing force exceeds the individual motor's stall force. Introduce the mapping (n+ , n− ) → n ≡ ( N+ + 1)n− + (n+ + 1) with 0 ≤ n ≤ N ≡ ( N+ + 1)( N− + 1). The corresponding probability density vector p satisfies equation (1.6) with the velocity in each state given by equation (1.8). The components Anm , n, m = 1, . . . , N of the state transition matrix A are given by the corresponding binding and unbinding rates of a given motor to a microtubule. Writing n ≡ n(n+ , n− ) the nonzero off-diagonal terms are [18, 89] Anm = π+ (n+ − 1), m = n(n+ − 1, n− ), (1.9a) Anm = π− (n− − 1), m = n ( n + , n − − 1), (1.9b) Anm = γ+ (n+ + 1), m = n(n+ + 1, n− ), (1.9c) Anm = γ− (n− + 1), m = n ( n + , n − + 1), (1.9d) π ( n ) ≡ ( N − n ) π0 , γ(n, F ) ≡ nγ0 e nFd , with F (1.10) where Fd is the experimentally measured force scale on which unbinding occurs. The diagonal terms are then Ann = − ∑m6=n Amn . We thus obtain a PDMP description of a tug-of-war model of cargo transport along microtubules. 1.2.1 Adiabatic Approximation In many applications, one finds that transitions between states are fast compared to v/∆, where v = maxn |vn | and ∆ is the fundamental length scale of the corresponding 12 problem. In our case, ∆ could represent the typical length of an axon, which ranges from 100-1000 microns. Performing the rescalings x → x/∆ and t → tv/∆ leads to the nondimensionalized version of equation (1.6): 1 ∂p = L(p) + Ap, ∂t ε (1.11) where 0 < ε 1. In this parameter regime, there are several transitions between motile states while the spatial position of the particle hardly changes at all. This suggests that the system rapidly converges to a quasi-steady state pss which is perturbed as x slowly changes. Performing this adiabatic approximation (also called the quasi-steady-state approximation) allows for easier analysis of the PDMP by recasting the full system as a single Fokker-Planck equation (sometimes called a Smoluchowski equation). The transition matrix A is assumed to be irreducible and conservative so that ψ ≡ (1, 1, . . . , 1)T is a left null vector of A. Let pss be an element of the right null space of A such that ψ T pss = 1. Let p ≡ ψ T p and w ≡ p − Cpss so that ψ T w = 0. Multiplying equation (1.11) on the left by ψ T gives ∂p = ψT L(Cpss + w). ∂t (1.12) On the other hand, substituting p = Cpss + w into equation (1.11) gives 1 ∂p ss ∂w p + = A(pss p) + L(pss p + w). ∂t ∂t ε (1.13) Combining equations (1.12) and (1.13) yields ∂w 1 = A(pss p + w) + (In − pss ψT )L(pss p + w), ∂t ε (1.14) where In is the n × n identity. Introduce the asymptotic expansion w = w0 + εw1 + ε2 w2 + . . . . (1.15) Substituting into equation (1.14) and collecting O(ε−1 ) terms yields Aw0 = 0. Since w is in the orthogonal complement of the left nullspace of A, it follows that w0 = 0. Now collecting O(1) terms gives Aw1 = −(In − pss ψ T )L(pss p). (1.16) Because ψ T (In − pss ψ T ) = 0, it follows that the right-hand side in the above equation is orthogonal to the null space of AT . The Fredholm Alternative Theorem guarantees the 13 existence of a solution for w1 , though it may not be unique. We obtain uniqueness by imposing ψ T w1 = 0. Hence, w ∼ εw1 . Substitution into equation (1.12) yields the Fokker- Planck (FP) equation : ∂p ∂ ∂ ∂p = − (V p ) + ε D , ∂t ∂x ∂x ∂x (1.17) V ( x ) = vT pss , (1.18) where the drift term V and diffusion coefficient D are given by D ( x ) = zT v, where v is the N-dimensional vector containing the velocity of each state and z an Ndimensional vector whose nth component is the solution to N ∑ m =1 Anm zm = [V ( x ) − vn ]pss n, (1.19) with ψ T z = 0. Equation (1.17) describes the time evolution of the probability density for the position of a particle evolving according to the PDMP in equation (1.6). If there is a sufficiently large number of particles whose probability densities evolve according to equation (1.17), then we may interpret equation (1.17) as a partial differential equation description of a bulk of motors moving in the domain. In the simple three-state model, we obtain an advection-diffusion equation for motion of bulk motors, with V= 1 v 1 − , Λ β+ β− D= ε (1 − V )2 (1 + V )2 + , Λ β2+ β2− Λ= 1 1 1 + + . β+ β− α (1.20) Hence, an advection-diffusion equation is our canonical mathematical representation of bulk motor dynamics. In this dissertation, we couple this equation with other processes to describe vesicular delivery and cellular length sensing. 1.3 Chemical Master Equation Throughout this dissertation, we will encounter situations where we have to model chemical reactions of N independent, identical structures with each structure's dynamics described by a continuous-time Markov process. In the limit N → ∞, the fraction of the population in a given state evolves according to a set of deterministic differential equations (also called kinetic equations). For finite N, one can track the stochastic fraction of structures in a given state using a chemical master equation. Let n = (n1 , n2 , . . . , nm ) denote the number of structures in each of m internal states with ∑m j=1 n j = N. The 14 probability that the population is in a state n at time t evolves according a master equation of the form dPn = dt ∑0 [Wnn0 Pn0 (t) − Wn0 n Pn (t)], (1.21) n with Wjk representing the reaction rate from the state j to state k. In general, equation (1.21) is very difficult to analyze on its own. However, it is possible to carry out a perturbation expansion of the master equation in terms of the system size 1/N and characterize the dynamics in terms of an equivalent Fokker-Planck equation or Langevin equation. Chemical master equations can be numerically simulated using the Gillespie algorithm [35]. 1.4 Structure of Dissertation In Chapter 2 and Chapter 3, we investigate motor-based vesicular delivery. Motivated by experimental observations in [78, 125], we mathematically investigate how allowing for reversible delivery of cargo to a particular target allows for a uniform distribution of vesicles. In Chapter 2, we investigate this in several geometries, including the semi-infinite interval, Cayley tree, disk, and sphere. In Chapter 3, we look at how including hard-core repulsion between individual motor-cargo complexes moving in bulk in our models impacts this phenomenon. In all the aforementioned cases, we see that irreversible delivery facilitates preferential vesicular delivery to targets proximal to the source of motors. Reversible delivery allows for uniformity in vesicular density provided motor velocities are not significantly hindered. In Chapter 4 and Chapter 5, we investigate how motor transport can be used for length control and length sensing in flagella and axons. In Chapter 4, we develop a stochastic version of a deterministic model for flagellar length control in [76], wherein the gradient of RanGTP at the base of a flagellum modulates the binding rate of intraflaggelar transport proteins (IFTs) at the base of the flagellum. We model the injection times of IFTs into the flagellum as a Poisson process whose rate in turn depends on a stochastic birth-death process describing the binding and unbinding of IFTs at the basal body. The result is a doubly stochastic Poisson process (DSPP) which models IFT injection into the flagellum. Furthermore, we invoke a Feynman-Kac formula to show that underlying the DSPP is a chemical master equation. In Chapter 5, we develop a mathematical version of a computational model for axonal 15 length sensing given in [102]. We show how axonal length is inversely related to the frequency of the oscillation of a somatic chemical signal carried from an axon's tip by dynein motors. We further show how the inverse relationship breaks down in the face of intracellular noise. We then describe a feed-forward network topology by which a cell may decode this information, and show how the decoding mechanism also breaks down in the face of noise. CHAPTER 2 CELL GEOMETRY AND VESICULAR DELIVERY A recent modeling study [17] investigated the active transport and delivery of vesicles across en passant synapses in the axons of neurons. Axons and dendrites both contain protein-rich synaptic subcellular compartments that form synaptic contacts between neurons. Some of the synaptic junctions occur at the terminals of axonal branches while others occur along the body of the axon. The latter are referred to as en passant synapses. The generation of new synaptic contacts during synaptogenesis or modification of old synapses in response to synaptic activity require localized protein delivery to a particular synaptic site [17, 61]. The relatively long distance between the soma and the distal axonal or dendritic synapses necessitates the use of active transport as a means of vesicular delivery. Hence, microtubules and molecular motors such as kinesin and dynein play key roles in proper synaptogenesis (see section 1.1). One main issue explored in reference [17] was the neuron's ability to evenly distribute vesicles across its en passant synapses-so-called synaptic democracy . Considering the fact that the source of motors that deliver cargo to these synapses is the soma of the neuron, one would expect that synapses proximal to the soma would obtain a greater amount of cargo compared to the distal axonal or dendritic synapses. However, the following experimental observations in axons of C. elegans and Drosophila [78, 125] suggest otherwise: (i) motordriven cargo exhibits ballistic anterograde or retrograde motion interspersed with periods of long pauses at presynaptic sites; (ii) the capture of vesicles by synapses during the pauses is reversible, in that vesicular aggregation at a site could be inhibited by signaling molecules resulting in dissociation from the target; (iii) the distribution of resources across synapses is relatively uniform. In reference [17], the transport and delivery of vesicles to synaptic targets was modeled using a one-dimensional (1D) advection-diffusion equation 17 based on the observations of [78, 125]. It was shown that in the case of irreversible cargo delivery, the steady-state vesicle density decays exponentially from the soma, whereas the steady-state density is relatively uniform in the reversible case. This suggests that reversibility in vesicular delivery plays a crucial role in achieving a "fair" distribution of resources within a cell. In this chapter, we significantly extend the 1D model of synaptic democracy in order to take into account the effects of cell geometry on reversible vesicular transport. We begin by briefly recounting the 1D results found in [17]; see section 2.1. Additionally, we investigate the behavior of the steady-state density of vesicles when the velocity of cargo-carrying motors is significantly different from free motors, which was not considered in reference [17]. We then consider a natural extension of the 1D analysis, namely a branching network (section 2.2). A tree is an appropriate domain to study synaptic democracy because it can account for the branched structure that is characteristic of axons and dendrites [88]. We show that in the irreversible case, branching increases the rate of decay of the steady-state distribution of vesicles. On the other hand, the steady-state profiles in the reversible case are similar to the 1D case. Moving away from highly polarized cells such as neurons, most cells (including a neuron's soma) have an approximately three-dimensional (3D) spherical shape. There are also examples of cells being treated as two-dimensional (2D) disks, particularly in the case of motile eukaryotic cells such as keratocytes [62, 80]. Therefore, we consider models of reversible vesicular transport in the disk and the sphere. We take the source of the motor-cargo complexes to be at the origin, and model the dynamics of the motor densities by differential equations transformed into their polar (2D) and spherical (3D) representations. The study of vesicular delivery in the disk and sphere domains is important because delivery of cargo to localized subcellular compartments is of utmost importance for several processes that occur in all cells. Such delivery is necessary, for example, when there is a need for restructuring a cell's cytoskeleton during cell growth, motility (see section 1.1.1), mitosis, and polarization [51] or when cellular waste materials are carried by autophagosomes to lysosomes for degradation [109]. In contrast to the 1D model, we distinguish between two types of filament distributions: (i) the distribution of microtubules emanating from the origin forms a continuum (section 2.3); (ii) the set of microtubules emanating from the origin forms a discrete set (section 2.4). In case (i), we 18 model the motion of motor densities using advection-diffusion equations. We find that for irreversible delivery, the steady-state vesicle density decays according to a modified Bessel function, whereas a uniform density can be obtained when delivery is reversible. In case (ii), we derive PDEs for the motor density based on stochastic differential equations (SDEs) for individual motor dynamics in the 2D and 3D domains following along the lines of Lawley et al. [72]. Throughout the chapter, we ignore boundary effects away from the source of motor-cargo complexes. In the case of exponentially decaying steady-state densities, this is a reasonable approximation provided that the spatial rate of decay is smaller than the size of the physical domain. 2.1 Semi-infinite Track Before elucidating our model and results, we briefly present the 1D results found in [17]. 2.1.1 Irreversible Delivery Consider a population of motor-cargo complexes or particles moving on a semi-infinite track, each of which carries a single synaptic vesicle precursor (SVP) to be delivered to a synaptic site. Assume that these particles are injected at the soma (x = 0) at a fixed rate J1 and that the distribution of synaptic sites along the axon is uniform. That is, at any given spatial point x, a particle can deliver its cargo to a synapse at a rate k. Neglecting interactions between particles, the dynamics of the motor-cargo complexes can be captured by the advection-diffusion equation [17] ∂u ∂2 u ∂u = −v + D 2 − ku, ∂t ∂x ∂x x ∈ (0, ∞), (2.1) where u( x, t) is the particle density along the microtubule track at position x at time t. Note that equation (2.1) can be derived from more detailed biophysical models of motor transport under the assumption that the rates at which motor-cargo complexes switch between different motile states are relatively fast [17, 89]; see section 1.2. In particular, the mean speed will depend on the relative times that the complex spends in different anterograde, stationary, and possibly retrograde states, whereas the diffusivity D reflects the underlying stochasticity of the motion. Equation (2.1) is supplemented by the boundary condition at x = 0: 19 J (u(0, t)) = J1 , J (u) ≡ vu − D ∂u . ∂x (2.2) Let c( x, t) denote the concentration of delivered vesicles to the presynaptic sites at x at time t with ∂c = ku − λc, ∂t (2.3) where λ denotes the degradation rate for vesicles. Note that in the irreversible delivery case, including vesicular degradation is necessary to prevent blowup in the solutions for c( x, t). This consideration is not necessary in the reversible delivery case. The steady-state solution for c is given by k J1 e−ξx , c= λ Dξ + v ξ≡ −v + √ v2 + 4Dk , 2D (2.4) which clearly indicates that c decays exponentially with respect to distance from the soma with correlation length ξ −1 . Taking the typical values D = 0.1µm s−1 for cytoplasmic diffusion and v = 0.1 − 1 µm s−1 for motor transport [48], and assuming that k 1 sec−1 , we see that ξ −1 ≈ (v/k) µm. Thus, in order to have correlation lengths comparable to axonal lengths of several millimeters, we would require delivery rates of the order k ∼ 10−5 sec−1 , whereas measured rates tend to be of the order of a few per minute [17, 45, 74]. This simple calculation establishes that injecting motor-complexes from the somatic end of the axon leads to an exponentially decaying distribution of synaptic resources along the axon. We now show, following reference [17], that relaxing the irreversible delivery condition in this model allows for a more uniform distribution of vesicles along the axon. 2.1.2 Reversible Delivery In order to take into account the reversibility of vesicular delivery to synapses, one must consider a generalization of the advection-diffusion model (2.1). To that end, let u0 ( x, t) and u1 ( x, t) denote the density of motor-cargo complexes without and with an attached SVP, respectively, and let k + and k − denote the rates at which vesicles are delivered to synaptic sites and recovered by the motors, respectively. Each density evolves according to an advection-diffusion equation combined with transition rates that represent the delivery and recovery of SVPs: 20 ∂u0 ∂u0 ∂2 u0 = − v0 + D 2 − γ0 u0 + k + u1 − k − cu0 , ∂t ∂x ∂x ∂u1 ∂u1 ∂2 u1 = − v1 + D 2 − γ1 u1 − k + u1 + k − cu0 , ∂t ∂x ∂x (2.5) (2.6) with x ∈ (0, ∞). Disparity in the velocities in each state reflects the effect cargo can have on particle motility, whilst the degradation rates γ0,1 are included to account for the possibility of particle degradation or recycling. Equations (2.5) and (2.6) are supplemented by the boundary conditions J (u j (0, t)) = Jj , j = 0, 1, (2.7) where Jj is the constant rate at which particles with or without cargo are injected into the axon from the soma. The dynamics for c( x, t) are now given by ∂c = k + u1 − k − cu0 . ∂t (2.8) We need not explicitly include degradation in this case because, provided J0 > 0, c( x, t) will be bounded. The steady-state distribution of vesicles is then c= k + u1 . k − u0 Substitution into the steady-state analogs of equations (2.5) and (2.6) yields q −ξ j − v + v2j + 4Dγ j j Jj e x , ξj = , u j (x) = Dξ j + v j 2D (2.9) (2.10) whence c= k + J1 Dξ 0 + v0 −Γx e , k − J0 Dξ 1 + v1 (2.11) with Γ ≡ ξ 1 − ξ 0 . It is evident that if Γ = 0, then c has a spatially uniform distribution. Suppose that the diffusion and degradation rates of motors do not change when carrying cargo. Then Γ = 0 would imply that the velocities of the cargo-carrying motors are equal to the velocities of the free motors. However, we would expect v1 < v0 due to the added load of the cargo on the motor, and that this would lead to a loss of synaptic democracy since Γ > 0. Indeed, values of v1 less than v0 lead to steady-state profiles of vesicle density reminiscent of the exponential decay behavior of the irreversible delivery case (see Figure 2.1), although the spatial rate of decay is mitigated by the presence of reversible delivery. Hence, attaining synaptic democracy also depends on physical properties of the cargo being carried. Large cargo, for example, may not be uniformly distributed throughout an axon whereas smaller cargo will. 21 1 v = v0 1 v = 0.75v0 1 v1 = 0.5v0 v1 = 0.25v0 v1 = 0.1v Vesicle Density 0.8 0.6 0 0.4 0.2 0 0 20 40 60 80 100 x [μm] Figure 2.1. Figure depicting the loss of synaptic democracy as disparity in velocities between free motors and cargo-carrying motors grows normalized so all curves fit in one frame. Parameter values are D = 0.1µm2 s−1 , γ0,1 = 0.01s−1 , k ± = 0.01s−1 , J0 = J1 , v0 = 0.1µms−1 . Vesicle density is normalized so that c(0) = 1. 2.2 Cayley Tree One limitation of the above model is that it does not capture the highly branched nature of an axon. Therefore, we now investigate irreversible and reversible delivery of vesicles to synapses on a tree. For simplicity, we consider an unbounded, regular tree Λ radiating from a unique origin with branching number z and segment length L (a Cayley tree); see Figure 2.2. We denote the origin, or the mother node, by α and the tree node opposite of the mother node by β. Let S1 be the set of z nodes connected to β. Similarly, let S2 consist of the z2 nodes that are connected to the vertices of the first generation and so on. The nth generation thus consists of zn nodes. Since all nodes (and their associated branches) of a given generation are equivalent for a regular tree, we can consider a single direct path through the tree and label the branch linking the node in Si−1 to the node in Si by i, i = 0, 1, 2, . . ., where S0 , S−1 are identified with the nodes β and α, respectively. Consider a population of motor-cargo complexes or particles moving on Λ, each of which carries a single synaptic vesicle precursor (SVP) to be delivered to a synaptic site. Motors are injected into the tree at a constant rate J1 at the mother node, α. Each branch is of finite length L, and we denote the point on each branch closest to α as x = 0 and 22 S1 S2 S0 α β Figure 2.2. Cayley tree Λ with z = 2. the point farthest away from α by x = L. The movement of the motors along a branch preceding a node in Si can be modeled by an advection-diffusion equation ∂ui ∂u ∂2 u = −v i + D 2i , ∂t ∂x ∂x (2.12) where ui ( x, t) represents the motor density at position x at time t, D is the motor diffusion coefficient, and v is the motor velocity. In the following, equation (2.12) will be coupled with the boundary conditions ui ( L, t) = ui+1 (0, t), i ≥ 0, J0 (0, t) = J1 , Ji ( L, t) = zJi+1 (0, t), i ≥ 1. (2.13) The first boundary condition represents continuity of motor density at the nodes of the tree. The second boundary condition represents the constant injection rate of motors at the mother node, and the last boundary condition reflects Kirchoff's law of conservation of current. Here, Ji ( x, t) = vui − D ∂ui . ∂x (2.14) Note that for simplicity, we take the motor velocity and diffusivities to be the same in all branches of the tree. A more detailed model would need to take into account a number 23 of features. For example, exclusion effects could mean motor velocities are locally density dependent, and diffusivities could change if the cross-sectional area of the axon decreases along the tree. Let us now use this setup to investigate irreversible and reversible vesicular delivery, respectively, to target synapses. 2.2.1 Irreversible Delivery We modify equation (2.12) by including a degradation term to account for irreversible delivery of vesicles. Let ci ( x, t) denote the concentration of vesicles at position x at time t on the i-th branch. The model for motor and vesicle dynamics is given by ∂ui ∂u ∂2 u = −v i + D 2i − kui , ∂t ∂x ∂x ∂ci = kui − λci , ∂t (2.15) (2.16) where λ is the vesicular degradation rate. At steady state, we have ∂ui ∂2 u + D 2i − kui = 0, ∂x ∂x kui c= . λ −v (2.17) (2.18) The general solution to equation (2.17) is given by ui ( x ) = Ai e ξ+ x + Bi e ξ− x , ξ± ≡ v± √ v2 + 4Dk , 2D (2.19) where Ai , Bi are constants of integration to be determined from boundary conditions. We can determine one of the constants for u0 by imposing the boundary condition reflecting the injection rate of motors. For the remaining constants, we employ the following method. We assume the motor density at each node in Si is given by Φi+1 . This ensures the solution on the tree will be continuous at the nodes. We then impose the boundary condition reflecting Kirchoff's law to determine each value Φi . That is, assume u0 (0) = Φ0 , (2.20) u0 ( L ) = u1 (0) = Φ1 , (2.21) u i −1 ( L ) = u i (0 ) = Φ i , i ≥ 2, (2.22) From equations (2.19) and (2.22), we have for i ≥ 1 ui ( x ) = Φ i e ξ − L − Φ i +1 ξ + x Φ i +1 − Φ i e ξ + L ξ − x e + ξ− L e . eξ − L − eξ + L e − eξ + L (2.23) 24 Imposing the current conservation condition (2.13), we obtain the following linear homogenous recurrence relation: z ( ξ + − ξ − ) Φ i +1 + ( ξ + − ξ − ) e ( ξ + + ξ − ) L Φ i −1 v(z − 1)(eξ − L − eξ + L ) + ξ − (eξ − L + zeξ + L ) − ξ + (eξ + L + zeξ − L ) Φi = 0. + D (2.24) Equation (2.24) has the solution Φi = νi with ν determined from the auxiliary equation z ( ξ + − ξ − ) ν2 + ( ξ + − ξ − ) e ( ξ + + ξ − ) L v(z − 1)(eξ − L − eξ + L ) + + ξ − (eξ − L + zeξ + L ) − ξ + (eξ + L + zeξ − L ) ν = 0. D (2.25) We obtain two solutions ν± in solving the quadratic equation, with |ν+ | > 1 and |ν− | < 1. Hence, i i Φi = c1 ν+ + c2 ν− . (2.26) In the case of an unbounded tree, we set c1 = 0, otherwise |Φn | → ∞ as n → ∞. Hence, i . Setting i = 0 gives c = Φ . Hence, we have Φi = cν− 0 i Φi = Φ0 ν− . (2.27) It remains to determine Φ0 . First, imposing the boundary conditions u0 ( L) = Φ1 and J0 (0) = J1 , the solution for u0 ( x ) is given by J1 eξ − L − (v − Dξ − )Φ1 eξ + x [v − Dξ + ]eξ − L − [v − Dξ + ]eξ − L u0 ( x ) = + (v − Dξ + )Φ1 − J1 eξ + L eξ − x . [v − Dξ + ]eξ − L − [v − Dξ − ]eξ + L (2.28) From equation (2.27), we obtain that Φ1 = Φ0 ν− . On the other hand, by substituting x = 0 into equation (2.28), we obtain Φ0 [v − Dξ + ]eξ − L − [v − Dξ − ]eξ + L − J1 (eξ − L − eξ + L ) Φ1 = . D (ξ − − ξ + ) (2.29) Equating the above two equations for Φ1 gives the explicit formula for Φ0 , Φ0 = − J1 eξ − L − eξ + L . Dν− (ξ − − ξ + ) − ([v − Dξ + ]eξ − L − [v − Dξ − ]eξ + L ) (2.30) We can now use equation (2.27) to obtain Φi , ∀i ∈ Λ. Hence, we have the steady-state distribution of vesicles in the Cayley tree in the irreversible delivery case. 25 In Figure 2.3, we compare the decay of vesicle density in the irreversible case of the Cayley tree to the semi-infinite track. We can see that toward the soma, the profiles are in exact agreement, whereas as soon as we reach the first branching point of the tree, the steady-state vesicle density suddenly drops, thereby aggravating the decay in the case of the Cayley tree to be greater than in the semi-infinite track. This suggests that if vesicular delivery were irreversible, biased delivery toward the soma would be greater than predicted in reference [17]. 2.2.2 Reversible Delivery To allow for re-uptake of vesicles from target sites, we must include the dynamics of cargo-carrying motors, u1i ( x, t) as well as free motors, u0i ( x, t), on each branch i and add switching terms to the advection-diffusion equation (2.12). Let ci ( x, t) represent the density of vesicles at position x at time t on branch i, i ∈ Sn . Then the motor and vesicle dynamics are given by ∂ui ∂2 u i ∂u0i = −v0 0 + D 20 − γ0 u0i + k + u1i − k − ci u0i , ∂t ∂x ∂x i i ∂u1 ∂u ∂2 u i = −v1 1 + D 21 − γ1 u1i − k + u1i + k − ci u0i , ∂t ∂x ∂x ∂ci = k + u1i − k − ci u0i . ∂t (2.31) (2.32) (2.33) We couple equations (2.31) and (2.32) with the boundary conditions (2.13). Let J0,1 be the injection rates at the origin for u0,1 , respectively. At steady state, we have ∂2 u0i ∂u0i − v − γ0 u0i = 0, 0 ∂x2 ∂x ∂2 u i ∂ui D 21 − v1 1 − γ1 u1i = 0, ∂x ∂x i k + u1 ci = . k − u0i D (2.34) (2.35) (2.36) As equations (2.34) and (2.35) are decoupled, we may apply the method elaborated in the irreversible delivery case to each equation separately and then obtain the full solution on the Cayley Tree. For example, let u00 (0) = Φ0 , u01 (0) = Ψ0 , (2.37) u00 ( L) = u10 (0) = Φ1 , u01 ( L) = u11 (0) = Ψ1 , (2.38) u0i−1 ( L) = u0i (0) = Φi , u1i−1 ( L) = u1i (0) = Ψi . (2.39) 26 10 Cayley Tree 9 1D Track Vesicle Density 8 7 6 5 4 3 2 1 0 0 20 40 60 80 100 x [μm] Figure 2.3. Plot comparing steady-state vesicle densities in the irreversible delivery case of the semi-infinite track and the Cayley tree. Parameter values are L = 10µm, v = 0.1µms−1 , D = 0.1µm2 s−1 , z = 3, λ = k = 0.01s−1 . Vesicle density is normalized so that c(0) = 1. so that the solutions are continuous at all nodes of Λ. We then have for i ≥ 1 Φ i e ξ − L − Φ i +1 ξ + x Φ i +1 − Φ i e ξ + L ξ − x + ξ− L e e , eξ − L − eξ + L e − eξ + L Ψ i e ζ − L − Ψ i +1 ζ + x Ψ i +1 − Ψ i e ζ + L ζ − x u1i ( x ) = ζ L + ζ− L e e , e − − eζ + L e − eζ + L u0i ( x ) = with ξ± ≡ v0 ± q v20 + 4Dγ0 2D , ζ± ≡ v1 ± q v21 + 4Dγ1 2D (2.40) (2.41) . (2.42) By imposing conservation of current at each node, we obtain second-order iterative equations for Φi and Ψi , which can be solved to give i Φi = ν− Φ0 , Ψi = µi− Ψ0 , (2.43) where ν− is the smaller root of equation (2.25) and µ− is the smaller root of the corresponding quadratic equation obtained by the replacement ξ ± → ζ ± . Finally, solving the equa- tions for u00 ( x ) and u01 ( x ) and imposing the boundary conditions u00 ( L) = Φ1 , u01 ( L) = Ψ1 yields Φ0 = − eξ − L − eξ + L J0 , D ν− (ξ − − ξ + ) − (ξ − eξ − L − ξ + eξ + L ) (2.44) 27 and Ψ0 = − J1 eζ − L − eζ + L . D ν− (ζ − − ζ + ) − (ζ − eζ − L − ζ + eζ + L ) (2.45) Thus, we obtain the steady-state vesicle distribution. Note that if γ0 = γ1 , then ξ ± = ζ ± when v0 = v1 and u1i ( x ) = u0i ( x ) for all x ∈ [0, L] and i ≥ 0. It follows that the vesicle distribution is uniform. In Figure 2.4, we show how uniformity in vesicle distribution is lost when v1 < v0 and compare decays in the Cayley tree with z = 3, L = 1 to decays on the semi-infinite track. We can see that the steady-state profiles are similar. Interestingly, the tree domain seems to facilitate a higher density of vesicles farther along in the domain than the semi-infinite track. We now investigate the impact of reversible vesicle delivery in higher-dimensional domains. 2.3 Higher-dimensional Geometries Although a 1D model is a reasonable first approximation of microtubule-based active transport in the axons and dendrites of a highly polarized cell such as a neuron, in most cells, intracellular transport takes place along 2D or 3D cytoskeletal networks of microtubules. For a sufficiently dense network, one could imagine carrying out some form of homogenization to obtain a continuum of microtubules. On the other hand, for a sparse network, the discrete nature of microtubules has to be taken into account. Here we focus on the continuum case; discrete microtubular networks will be considered in section 2.4. For simplicity, we model a cell as a disk or a sphere and assume that the density of microtubules is radially symmetric, that is, we ignore the curvature of microtubules. We take the source of the motor-cargo complexes to be at the origin of the cell, and represent the dynamics of the motor densities by advection-diffusion equations transformed into their polar (2D) and spherical (3D) representations. We will also assume that each motor carries one cargo element and can deliver its cargo at any point within the given domain. In other words, we assume that there is a continuum of target sites within the cell. 28 1 v =v 1 0 v = 0.75v 1 0 v1 = 0.5v0 v1 = 0.25v0 v1 = 0.1v0 Vesicle Density 0.8 0.6 0.4 0.2 0 0 20 40 60 80 100 x [μm] Figure 2.4. Plots showing loss of vesicular uniformity as v1 decreases in the case of a Cayley tree (kinked curves) and a semi-infinite track (smooth dotted curves). Parameter values are v0 = 0.1µs−1 , D = 0.1µ2 s−1 , γ0,1 = 0.01s−1 , k ± = 0.01s−1 , L = 10µm, z = 3, J0 = J1 . Vesicle density is normalized so that c(0) = 1. 2.3.1 The Disk Let Ω2 ≡ R2 \ Bδ (0), where Bδ (0) is the disk of radius δ centered at the origin, with 0 < δ 1. In polar coordinates, Ω2 = {(r, θ )|r ≥ δ, 0 ≤ θ ≤ 2π }. (2.46) We model the dynamics of the motor population by an advection-diffusion equation that is a radially-symmetric 2D analog of the 1D model. As in the previous cases, we first consider irreversible vesicle delivery and then reversible vesicle delivery. Irreversible Delivery. Let u(r, t) and c(r, t) denote, respectively, the density of motors and vesicles at a radial distance r from the origin at time t. The motor and vesicle densities are taken to evolve according to the equations ∂u ∂2 u D − V ∂u =D 2 + − ku, ∂t ∂r r ∂r ∂c = ku − λc, ∂t (2.47) (2.48) 29 where D is the diffusion coefficient, v = V /r is a divergence-free motor velocity∗ , and λ is the degradation rate of vesicles. As in the 1D case, we model irreversible vesicle delivery using an effective degradation term in equation (2.47). We pair equation (2.47) with the boundary conditions u ( δ ) = u0 , lim u(r ) = 0, r →∞ (2.49) where u0 > 0 denotes the density of motors on ∂Bδ (0). At steady state, we have the equations ∂2 u D − V ∂u + − ku = 0, ∂r2 r ∂r ku c= . λ D (2.50) (2.51) The steady-state vesicle density profile is immediately given by a modified Bessel function of the second kind: r 2D ku0 D/k c (r ) = δ . V λ 2D δ KV √ 2D D/k V r 2D K V √ (2.52) As in previous geometries, irreversible vesicular delivery results in a decaying steady-state profile for vesicle density. In Figure 2.5, we compare the decay in the disk with the decay on the semi-infinite track. We can see that towards the origin, the Bessel function distributes vesicles more liberally than the exponential function but then rapidly decays below the latter. We also show a plot of the corresponding decay in the case of a sphere (see section 2.3.2), which is similar to the disk. Let us now look at the reversible vesicle delivery case. Reversible Delivery. To account for the possibility of re-uptake of vesicles by free motors, we model the dynamics of both the free motor density, u0 (r, t), and the cargo-carrying motor density, u1 (r, t). We thus have a pair of radially-symmetric advection-diffusion equations coupled with switching terms that reflect vesicle delivery and uptake. Again, let c(r, t) ∗ This is motivated by the idea that the density of microtubules decreases as r −1 in the 2D case (and decreases as r −2 in the 3D case). When we consider a discrete distribution of microtubules, the effective velocity will have a more complicated r-dependence 30 10 1D 2D 3D 9 Vesicle Density 8 7 6 5 4 3 2 1 0 0 20 40 60 80 100 x [μm] Figure 2.5. Figure comparing irreversible vesicular profiles from equations (2.52) and (2.4). Here x either represents 1D distance or a radial coordinate. Parameter values are D = 0.1µ2 s−1 , V = 1µ2 s−1 , λ = k = 0.01s−1 , δ = 0.1µm, and the flux J1 chosen appropriately so as to match up to the left boundary data. Also shown is the corresponding steady-state density for the sphere. Vesicle density is normalized so that c(0) = 1. denote the vesicle density at a distance r from the origin at time t. The system of equations is ∂u0 ∂2 u0 D ∂u0 =D 2 + − ∂t ∂r r ∂r ∂2 u D ∂u1 ∂u1 = D 21 + − ∂t ∂r r ∂r ∂c = −k + cu0 + k − u1 , ∂t V0 ∂u0 − γ0 u0 − k + cu0 + k − u1 , r ∂r V1 ∂u1 − γ1 u1 + k + cu0 − k − u1 , r ∂r (2.53) (2.54) (2.55) where D is the motor diffusion coefficient, v0,1 = V0,1 /r are divergence-free velocities of the free and cargo-carrying motors, respectively, k ± denote the rates of vesicle uptake and delivery, respectively, and γ0,1 are motor degradation rates. We again point out that the reversibility in vesicle delivery means that we do not need to include a degradation term in equation (2.55). The corresponding system at steady state is ∂2 u0 D ∂u0 V0 ∂u0 + − − γ0 u0 = 0, 2 ∂r r ∂r r ∂r ∂2 u1 D ∂u1 V ∂u D 2 + − 1 1 − γ1 u1 = 0, ∂r r ∂r r ∂r k − u1 c= , k + u0 D (2.56) (2.57) (2.58) 31 The solution for ul (r, t), l = 0, 1, is r γ l r K Vl r D 2D 0 ul = ul r γ , Vl l D δ δ K Vl D 2D Vl 2D (2.59) where u0l denotes the boundary data at the origin for ul . It immediately follows that r γ r γ 0 1 V −V K V δ K V1 r k − u01 r 12D 0 2D0 D D 2D r r c= (2.60) γ γ . k + u00 δ 0 1 K V1 δ K V0 r D D 2D 2D Assume that motor degradation rates are equal, γ1 = γ0 . It is clear that if V1 = V0 , then the vesicle distribution is uniform. If V1 < V0 , then the spatial profile is a decaying function of r; see Figure 2.6(a). The behavior here is consistent with what is seen along the semi-infinite track, although in the latter case, the decay is exponential. As expected, the rate of decay is mitigated by a reduction in the motor degradation rates as shown in Figure 2.6(b). 2.3.2 The Sphere Let Ω3 ≡ R3 \ Bδ (0), where Bδ (0) is the ball of radius δ centered at the origin, with 0 < δ 1. In spherical coordinates, the domain is defined as Ω3 = {(ρ, θ, φ)|ρ ≥ δ, 0 ≤ θ ≤ 2π, 0 ≤ φ ≤ π }. (2.61) As in the case of a disk, we consider a population of motors sourced at the origin switching between diffusive and ballistic transport, depending on whether or not a given motor is bound to a microtubule. The dynamics of the motor population is modeled by a radiallysymmetric 3D advection-diffusion equation analogous to the 1D model. Let u(ρ, t) denote the density of motors located at a radial distance r from the origin at time t. Irreversible Delivery. Let c(ρ, t) represent the density of vesicles at a distance of ρ from the origin at time t. Then the dynamics of the motor and vesicle densities are given by ∂u ∂2 u 2D V ∂u =D 2+ − 2 − ku, (2.62) ∂t ∂ρ ρ ρ ∂ρ ∂c = ku − λc, (2.63) ∂t 32 Effects of cell geometry on reversible vesicular transport 14 1 Vesicle Density 0.8 V1 = V 0 V1 = 0.25V0 V1 = 0.75V0 V1 = 0.01V0 V1 = 0.50V0 0.6 0.4 0.2 (a) 0 0 20 40 60 80 100 x [μm] 1 Vesicle Density 0.8 0.6 0.4 γ 0.2 =10-3s-1 γ =2.5x10-3 s-1 0,1 0,1 γ =7.5x10-3 s-1 γ =10-2 s-1 0,1 γ =0 0,1 0,1 (b) 0 0 20 40 60 80 100 x [μm] Figure Figure 6: Figure showing loss of vesicular uniformity on the disk as V decreases. (a) 2.6. Figure showing loss of vesicular uniformity as V1 decreases.1 (a) Plot of Plot of steady-state steady-state vesicle density for various V1 values andmotor fixeddegradation motor degradation vesicle density for various V1 values and fixed rates −1 − 1 2 s−1 . and rates γ0,1 =0.01s 0.01s. (b). Corresponding (b) Corresponding plots degradation for various degradation rates γ0,1 = plots for various rates and V1 = 0.75µ 2 s−1 ,D = 0.1µm2 s−1 , δ =2 0.1µm, 2 −1 Other parameter values are V0 = 1µm 0.01s−12, s−1 , V1 = 0.75µm s . Other parameter values are V0 = 1µm s−1 , Dk± == 0.1µm 0 0 u0 = uk1 . = 0.011s−1 , γ −1 δ = 0.1µm, , u00 = u01 . ± 0,1 = 0.10s if V1 < V0 , then the spatial profile is a decaying function of r, see Fig. 6(a). The behavior here is consistent with what is seen along the semi-infinite track, although in the latter case the decay is exponential. As expected the rate of decay is mitigated by a reduction in the motor degradation rates as shown in Fig. 6(b). 33 where D is the motor diffusion coefficient, V /ρ2 is a divergence-free motor velocity, and λ is the vesicular degradation rate. As in the previous analysis, vesicular degradation must be accounted for in the irreversible delivery case to ensure vesicle profiles do not blow up. It is not necessary in the reversible case. We pair equation (2.62) with the boundary conditions u ( δ ) = u0 , lim u(ρ) = 0, ρ→∞ where u0 > 0 is the density of motors on ∂B. At steady state, we have the system 2D V ∂u ∂2 u − 2 − ku = 0, D 2+ ∂ρ ρ ρ ∂ρ ku c= . λ (2.64) (2.65) (2.66) As the steady-state equations are difficult to solve analytically, we solve them numerically. In Figure 2.5, we compare the decay of the steady-state vesicle density in 3D with the 1D and 2D domains. We find that the 3D steady-state profile behaves similarly to the 2D case. Reversible Delivery. In the reversible delivery case, we keep track of the free motor densities, u0 (ρ, t) and the cargo-carrying motor densities, u1 (ρ, t). We model the motor dynamics with advection diffusion equations coupled with switching terms to reflect delivery and uptake of vesicles to and from target sites. Let c(ρ, t) denote the vesicle density at a distance ρ from the origin at time t. The system capturing the dynamics is D ∂ 2 ∂u0 ∂u0 = 2 ρ − ∂t ρ ∂ρ ∂ρ ∂u1 D ∂ 2 ∂u1 = 2 ρ − ∂t ρ ∂ρ ∂ρ ∂c = k − u1 − k + cu0 , ∂t V0 ∂u0 − γ0 u0 − k + cu0 + k − u1 , ρ ∂ρ V1 ∂u1 − γ1 u1 + k + cu0 − k − u1 , ρ ∂ρ (2.67) (2.68) (2.69) where the various parameters are as in previous examples.. At steady state, we have the system, D ∂2 u0 + ∂ρ2 2D V0 − 2 ρ ρ ∂u0 − γ0 u0 = 0, ∂ρ ∂2 u1 2D V1 ∂u1 D 2 + − 2 − γ1 u1 = 0, ∂ρ ρ ρ ∂ρ k − u1 c= . k + u0 (2.70) (2.71) (2.72) 34 Again, we obtain the steady-state profiles numerically. Clearly, if V1 = V0 , we have a uniform distribution of vesicles. When V1 < V0 , we again have similar behavior to the 2D profiles; see Figure 2.7. An explicit comparison of the distributions in the 1D, 2D, and 3D cases is shown in Figure 2.8. 2.4 Discrete Microtubule Distributions The models in section 2.3 were phenomenologically based, under the assumption that we could treat a cytoskeletal network as a continuum, and model the effective motor transport as a radially-symmetric advection-diffusion equation. It is possible to extend the adiabatic analysis of section 1.2.1 in order to derive a higher-dimensional advection-diffusion equation from a more realistic stochastic model of 2D or 3D motor transport, in which individual motors switch between ballistic motion when bound to a microtubule and diffusive motion when unbound [19]. In general, the resulting advection-diffusion equation will be anisotropic, with an associated diffusion tensor that depends on the configuration of microtubules. Here we will consider a different regime in which the cytoskeletal network is sparse so that we have a discrete network. In order to simplify our analysis, we will assume that the microtubules project radially from the center of the disk or sphere. We can then derive an effective advection-diffusion equation for motor transport by following recent analysis of virus trafficking in cells [69, 72]. 2.4.1 The Disk Consider a finite set of N identical, evenly spaced microtubules radiating from the center of the disk [69, 72]. That is, Ω2 is partitioned into N equal slices, each of angular width Υ ≡ 2π/N (see Figure 2.9), whose boundaries correspond to microtubules. Following Lawley et al. [72], we will derive an effective advection-diffusion equation for motor transport by considering the dynamics of a single molecular motor moving within a single slice U2 ≡ [δ, ∞) × [0, Υ] ⊂ Ω2 - restriction to a single slice is allowed because of the symmetric partitioning and the fact that we are only interested in the radial distribution of motors. Therefore, consider a single motor-cargo complex originating on ∂Bδ and undergoing Brownian motion in the interior of U2 until it reaches a microtubule, whence it binds 35 Effects of cell geometry on reversible vesicular transport 16 1 Vesicle Density 0.8 V1 = V0 V1 = 0.25V0 V1 = 0.75V0 V1 = 0.10V0 V1 = 0.50V0 0.6 0.4 0.2 (a) 0 0 20 40 60 80 100 x [μm] 1 Vesicle Density 0.8 0.6 0.4 γ 0,1 =7.5x10-3 s-1 γ =10-2 s-1 0,1 γ =0 =10-3s-1 γ =2.5x10-3 s-1 0,1 γ 0.2 0,1 0,1 (b) 0 0 20 40 60 80 100 x [μm] Figure 7: Figure depicting loss of uniform vesicle distribution on the sphere when V1 Figure 2.7. Figure depicting loss of uniform vesicle distribution on the sphere when V1 decreases.decreases. (a) Steady state distributions of vesicles inin3D 3Dsphere sphere for various V1 values (a) Steady-state distributions of vesicles for various V1 values −1 and fixedand motor ratesrates γ0,1γ0,1 == 0.01s Corresponding for various fixeddegradation motor degradation 0.01s−1.. (b) (b) Corresponding plotsplots for various 3 3 −1 degradation ratesVand = 0.75µmss−1 .. Other parameter valuesvalues are as inare Figure degradation rates and Other parameter as 2.6. in Fig. 6. 1 =V10.75µm ! " 2D V1 ∂u1 ∂ 2 u1 + − 2 − γ 1 u1 = 0 D ∂ρ2 ρ ρ ∂ρ k− u1 c= . k + u0 Effects of cell geometry on reversible vesicular transport 17 36 1 0.8 Vesicle Density 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 0.6 0.4 0.2 0 0 20 40 60 80 100 x [μm] Figure 8: Comparison of profiles in 1D (dotted), 2D (dashed), and 3D (solid) domains −1 −1 (dotted), 2D (dashed), and 3D (solid) Figure Comparison profiles in2 s1D domains for v1 2.8. = 0.075µms , Vof 0.75µm for the disk and V1 = 0.75µm3 s−1 for the 1 = − 1 2 − 1 3 − 1 for v1 = Other 0.075µms , V1 =values 0.75µm the disk and 1 and 6.V1 = 0.75µm s for the sphere. sphere. parameter ares as for in Figs. Other parameter values are as in Figure 2.1 and Figure 2.6 5. Discrete Microtubule Distributions to the microtubule and moves ballistically away from the origin for some exponentially The models in section 4 were phenomenologically-based, under the assumption that distributed amount of time. At this point, the motor-cargo complex is reinserted into the we could treat a cytoskeletal network as a continuum, and model the effective motor transport a radially symmetric advection-diffusion equation. It is possible slice at theas current radius for some randomly selected angle between 0 andtoΥ.derive If X (t) a higher-dimensional advection-diffusion equation from a more realistic stochastic represents radialtransport, distance from the origin and θmotors (t) represents angle model of the 2D motor's or 3D motor in which individual switch some between ballistic [0, motion bound to aismicrotubule unbound between Υ], thewhen motor's motion described byand the diffusive followingmotion systemwhen of SDEs [69, 72]: [3]. In general, the resulting advection-diffusion equation will be anisotropic, with an associated diffusion tensor that depends on the configuration of microtubules. Here Vdt, θ = 0, Υ, √ dX = regime in which we will consider a different the cytoskeletal network is sparse so that ( D/X )dt + 2DdWX , θ ∈ (0, Υ), In order to simplify our analysis, we will assume that we have a discrete network. 0,√ θ = 0, Υ, the microtubules project from the center of the disk or sphere. We can then (2.73) dθ =radially ( 2D/X )dWθ , θfor∈ motor (0, Υ),transport by following recent derive an effective advection-diffusion equation analysis of virus trafficking in cells [12, 18]. where WX , Wθ are standard independent Wiener processes, V is the motor velocity, and Themotor disk diffusion coefficient. Note that one major difference from models of virus D5.1. is the Consider ais finite setare of N identical, spacedtransport microtubules radiating the at trafficking that we interested in evenly the outward of motors fromfrom a source center of the disk [12, 18]. That is, Ω2 is partitioned into N equal slices, each of the origin,width whereas enterFig. the 9), cellwhose at some finite distance R fromtothe cell center and angular Υ ≡viruses 2π/N (see boundaries correspond microtubules. we cell willnucleus. derive an equation Following Lawley et al. [18], the move inwards in order to find Ineffective referenceadvection-diffusion [72], Lawley et al. use a coarse for motor transport by considering the dynamics of a single molecular motor moving graining to derive a single effective describing the overall radial motion of within amethod single slice U2 ≡ [δ, ∞) × [0, Υ] ⊂ ΩSDE 2 - restriction to a single slice is allowed the symmetric partitioning and(2.73). the fact that we arethere onlyisinterested in the abecause particle of evolving according to equations They assume a continuous-time radial distribution of motors. jump Therefore, Markov process underlying themotor-cargo particle's switching diffusive andδ ballistic consider a single complexbetween originating on ∂B and undergoing Brownian motion in the interior of U2 until it reaches a microtubule, dynamics, and that the dynamics of the Markov process are very fast relative to all other whence it binds to the microtubule and moves ballistically away from the origin 37 Υ δ microtubules N=5 Figure 2.9. Partitioning of domain Ω2 for N = 5. processes. Invoking an adiabatic approximation, they derive the following coarse-grained effective SDE approximation to equations (2.73): dX = T(X) µ +V dt + X µ + T(X) µ + T(X) D s 2D T(X) dW, µ + T(X) (2.74) where W (t) is a standard Wiener process, µ is the mean for the exponential distribution dictating the amount of time a particle spends in the ballistic phase, and T ( X ) is the mean first passage time (MFPT) for a particle in the cytoplasm to reach a microtubule, T(X) = Υ2 X 2 . 12D (2.75) Let p(r, t) represent the probability that a particle evolving according to equation (2.74) is at a distance r from the origin at time t. The corresponding Fokker-Planck equation is i ∂p ∂ h D T (r ) µ ∂2 T (r ) =− +V p + 2 D p . ∂t ∂r r µ + T (r ) µ + T (r ) ∂r µ + T (r ) (2.76) Now suppose that there are N independent motors evolving according to the SDE (2.74). Let u(r, t) denote the density of motors at time t located a radial distance of r from the origin. We have the following relationship between p(r, t) and u(r, t): p(r, t) = where N = 2π 2π ru(r, t), N Z ∞ δ ru(r, t)dr. (2.77) (2.78) 38 We are assuming that u decays sufficiently fast at infinity. Substituting equation (2.77) into (2.76) yields the following PDE for motor density dynamics: i 1 ∂2 ∂u T (r ) 1 ∂ h T (r ) µ D =− D + Vr u + ru . ∂t r ∂r µ + T (r ) µ + T (r ) r ∂r2 µ + T (r ) (2.79) Using equation (2.79) as a starting point, we now investigate reversible vesicle delivery for the discrete microtubule set case. Consider the dynamics of free motors with density u0 (r, t) and cargo-carrying motors with density u1 (r, t). Each evolves according to an equation of the form (2.79), coupled with switching terms that reflect vesicle delivery and uptake. Again, let c(r, t) denote the vesicle density at a distance r from the origin at time t. The system of equations is then i 1 ∂2 T (r ) T (r ) 1 ∂ h µ ∂u0 D u0 + D =− + V0 r ru0 ∂t r ∂r µ + T (r ) µ + T (r ) r ∂r2 µ + T (r ) − k + cu0 + k − u1 , i 1 ∂2 ∂u1 T (r ) 1 ∂ h T (r ) µ D =− D + V1 r u1 + ru 1 ∂t r ∂r µ + T (r ) µ + T (r ) r ∂r2 µ + T (r ) + k + cu0 − k − u1 , ∂c = k − u1 − k + cu0 , ∂t (2.80) (2.81) (2.82) with D, V0,1 , k ± defined as in previous cases. One important difference is that we no longer include motor degradation terms, since these would lead to a breakdown of the analysis of Lawley et al. [72]. Again we find that a uniform steady-state distribution of vesicles occurs when V0 = V1 , but there is a loss in uniformity when V1 < V0 . This is illustrated in Figure 2.10(a) where we show numerical plots of the steady-state solutions. 2.4.2 The Sphere Let Ω3 be defined as in section 2.3.2, but now we let the set of microtubules emanating from the origin be discrete rather than a continuum; see Figure 2.11. Hence, we have a natural partition for Ω3 = ω1 ∪ ω2 . We define ω1 in the following way. Let N be the number of microtubules emanating out from the small sphere of radius δ enveloping the origin. These can be modeled as infinite cylinders each of radius ε. Let ci for i = 1...N denote a randomly selected fixed position on the δ-sphere. Then each microtubule Mi is defined as follows: n o Mi ≡ x ∈ Ω3 ||x − ρci || ≤ ε, ρ ∈ [δ, ∞) . (2.83) 39 Effects of cell geometry on reversible vesicular transport 20 1 Vesicle Density 0.8 0.6 0.4 0.2 0 (a) 0 20 40 60 80 100 x [μm] Vesicle Density 1 V1 = V0 V1 = 0.75V0 V1 = 0.50V0 0.8 V1 = 0.25V0 V1 = 0.10V0 0.6 0.4 0.2 0 0 (b) 20 40 60 80 100 x [μm] Figure 10: Figure showing the loss of vesicle uniformity for V1 < V0 = 1 and discrete Figure 2.10. Figure showing the loss of vesicle uniformity for V1 < V0 = 1 (b) and The discrete distribution of microtubules. (a) The disk with N = 12 microtubules. sphere distribution of microtubules. (a) The disk with N = 12 microtubules. 2(b)−1 The sphere with N = with 1000Nmicrotubules. OtherOther parameter values 0.1µm , δ0.1µm, = 0.1µm, 2 s−1s, δ = = 1000 microtubules. parameter valuesare are D D = = 0.1µm −1 0 −1 00 k± = 0.01sk± =, 0.01s u0 = u , u10. = u01 . number of microtubules emanating out from the small sphere of radius δ enveloping the origin. These can be modeled as infinite cylinders each of radius ε. Let ci for i = 1...N denote a randomly selected fixed position on the δ-sphere. Then each microtubule Mi is defined as follows: " ! # " 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 T (X) = X D − 2 ε ln N X + ln 4 − 1 − 4 Ψ N2 with Ψ= N N & & k=1 j=k+1 ln ||ck − cj || 40 δ Figure 11: Sketch of Ω3 showing N = 6 microtubules radiating from center. Figure 2.11. Sketch of Ω3 with N = 6. Let p(ρ, t) represent the probability that a particle is at position ρ at time t. The Fokker-Planck equation associated with equation (5.8) is We take ω1 = ∪iN= M and ω = Ω3 \ ω1 . To the dynamics "of motor-cargo com" model % ∂ $!1 2Di T (ρ) 2 µ ∂ 2 ! T (ρ) ∂p =− +V p +D 2 p (5.9) plexes in we the motion of a ∂t this domain, ∂ρ ρ µ + must T (ρ) derive µ + PDEs T (ρ) from the ∂ρ SDEs µ + describing T (ρ) Now suppose that there are N independent motors evolving according to equation single particle in this domain. WeS assume a single particle's motion is characterized by (5.8). If u(ρ, t) denotes motor density at a distance ρ from the origin at time t, we have Brownian the following relation t) anda microtubule, u(ρ, t) standard motion in ωbetween reaches when it undergoes ballistic 2 until itp(ρ, 4π 2 ρ from u(ρ, t)the origin for some exponentially distributed (5.10)time. p(ρ, t) V = away motion with fixed velocity NS where is then released at a random position in Ω3 with radius equal to how far it The particle ' ∞ ρ2 u(ρ,et t)dρ S = 4π reached with ballisticNmotion. Lawley al. provide the following SDE as a coarse-grain δ approximation to a particle moving through Ω3 as described above: s 2D T ( X ) T(X) µ dX = +V dt + 2D dW, X µ + T(X) µ + T(X) µ + T(X) (2.84) where X (t) is the distance of a particle from the origin, D is the diffusion coefficient, V is the particle velocity, W (t) is a Wiener process, µ is the mean for the exponential density for the amount of time a particle spends on a microtubule, and T ( X ) is the MFPT for a particle beginning at position X to reach a microtubule. Coombs, Straube, and Ward provide the following asymptotic approximation for T ( X ) in the small ε limit: T(X) = ε X2 h 2 4 i − ln + ln 4 − 1 − 2 Ψ , D N X N (2.85) with Ψ= N N ∑ ∑ k =1 j = k +1 ln ||ck − c j ||. (2.86) 41 Let p(ρ, t) represent the probability that a particle is at position ρ at time t. The Fokker -Planck equation associated with equation (2.84) is i ∂p ∂ h 2D T (ρ) µ ∂2 T ( ρ ) =− +V p +D 2 p . ∂t ∂ρ ρ µ + T (ρ) µ + T (ρ) ∂ρ µ + T (ρ) (2.87) Now suppose that there are NS independent motors evolving according to equation (2.84). If u(ρ, t) denotes motor density at a distance ρ from the origin at time t, we have the following relation between p(ρ, t) and u(ρ, t) p(ρ, t) = where NS = 4π 4π 2 ρ u(ρ, t), NS Z ∞ ρ2 u(ρ, t)dρ. (2.88) (2.89) δ Substituting equation (2.88) into (2.87) gives i D ∂2 T ( ρ ) ∂u 1 ∂ h T (ρ) µ =− 2 2Dρ + Vρ2 u + 2 2 ρ2 u . ∂t ρ ∂ρ µ + T (ρ) µ + T (ρ) ρ ∂ρ µ + T (ρ) (2.90) Equation (2.90) is the PDE describing the dynamics of motor density in Ω3 . We now use it as the governing PDE to investigate reversible vesicular delivery in a sphere. Consider the dynamics of free motors with density u0 (r, t) and cargo-carrying motors with density u1 (r, t). Each evolves according to an equation of the form (2.79), coupled with switching terms that reflect vesicle delivery and uptake. Again, let c(r, t) denote the vesicle density at a distance r from the origin at time t. The system of equations is then i D ∂2 T ( ρ ) ∂u0 1 ∂ h T (ρ) µ =− 2 2Dρ + V0 ρ2 u0 + 2 2 ρ2 u0 ∂t ρ ∂ρ µ + T (ρ) µ + T (ρ) ρ ∂ρ µ + T (ρ) − k − cu0 + k + u1 , (2.91) h i 2 ∂u1 1 ∂ T (ρ) µ D ∂ T (ρ) =− 2 2Dρ + V1 ρ2 u1 + 2 2 ρ2 u1 ∂t ρ ∂ρ µ + T (ρ) µ + T (ρ) ρ ∂ρ µ + T (ρ) + k − cu0 − k + u1 , (2.92) ∂c = k + u1 − k − cu0 . (2.93) ∂t Again we find that a uniform steady-state distribution of vesicles occurs when V0 = V1 , but there is a loss in uniformity when V1 < V0 . This is illustrated in Figure 2.10(b) where we show numerical plots of the steady-state solutions. 42 2.5 Discussion In this chapter, we investigated a possible biophysical mechanism for facilitating a more uniform distribution of vesicles to targets dispersed throughout cells of various geometries. In particular, we generalized the results found in [17] by examining the impact of allowing for reversibility in vesicular delivery to target sites on a Cayley tree, a disk, and a sphere. On the disk and sphere, we considered both a continuous distribution of microtubules and a discrete set. In the latter case, we derived an effective advection-diffusion equation for motor transport based on SDEs for single motor motion, following along similar lines to reference [72]. In all cases, we found that uniformity in the steady-state vesicle distribution is attainable if vesicle delivery is reversible, and the velocity of cargo-carrying motors is not significantly less than that of free motors. We also characterized the loss of uniformity when there was a mismatch between the velocities of free and cargo-bound motors. There are a number of possible extensions of this work (beyond taking into account exclusion effects as detailed in Chapter 3). For example, we assumed each motor carried only one cargo element. It would be interesting to relax this condition and allow each motor to carry a collection of cargo elements that can be delivered to target sites. This problem was previously investigated for a 1D track using aggregation theory and a modified version of the Becker-Doring equations in [13]. A related problem is developing a more detailed model of bidirectional motor transport. This is particularly important in determining how the effective speed of a motor-cargo complex depends on the cargo load, for in order to a achieve a more uniform distribution of resources using the proposed mechanism, it is necessary that the speed be weakly dependent on cargo load. In the case of large vesicles, this would require that transport involves cooperation between multiple molecular motors, rather than a single motor. There is considerable debate in the literature regarding the most likely mechanism for cooperative bidirectional transport [42, 68, 87]: (a) an asymmetric tug-of war model involving the joint action of multiple kinesin and dynein motors pulling in opposite directions; (b) a symmetric tug-of-war model where all the motors are of the same type, but they are distributed on microtubules of opposite polarity; (c) a hopping model, in which the whole motor-cargo complex hops between microtubules of opposite polarity; (d) some form of coordination complex that controls 43 the switching between different motor species. Yet another related issue is developing a more detailed model of the exchange of vesicular cargo between motors and targets. In this chapter, we simply took the exchange to be given by first-order kinetics, and assumed that there was a uniform distribution of targets throughout the cell. The latter is likely to be a particularly crude approximation in the case of higher-dimensional cell geometries. For example, in the case of discrete microtubular networks, one might expect targets to be located within some local neighborhood of the microtubules. A natural extension of the analysis on a tree would be to take successive generations of the tree to have different properties, reflecting the fact that in axons (and dendrites) the branches tend to taper off (become thinner). Another interesting problem is how one would extend the analysis of Lawley et al. [72] to more general configurations of microtubules. One of the essential steps in their analysis is the adiabatic approximation that during the time between binding and unbinding to a microtubule, the relative change in radial position is small. This has several implications for our own analysis. First, the adiabatic approximation breaks down at sufficiently large radii, as can be seen from the formula for the MFPT in equation (2.75), that is, T ( X ) ∼ X 2 . Thus, a more careful analysis would need to restrict the dynamics to a bounded domain and take the number of microtubules to be sufficiently large. The adiabatic requirement also meant that we had to neglect the degradation of motor-cargo complexes. Again, it would be interesting to extend the analysis of reference [72] to allow for the possibility that motors disappear so that one has to determine a conditional MFPT. CHAPTER 3 VESICLE DELIVERY WITH EXCLUSION In Chapter 2, we investigated the impact reversibility in vesicular delivery to target sites had on vesicular distribution throughout a cell's body. In particular, we showed that allowing for reversibility in vesicular delivery to target sites allows for uniform vesicular distribution provided the cargo-carrying motor's velocity is not significantly hampered by the load of the cargo. However, one inherent assumption in all the models discussed in Chapter 2 was that there was no interaction between individual motor-cargo complexes. In particular, one inherent assumption was that two motor-cargo complexes could occupy the same location at the same time. In this chapter, we generalize the results of reference [17] by considering the effects of exclusion between particles on the steady-state distribution of synaptic vesicles. As in [17], we compare the effects of reversible and irreversible vesicular delivery on this distribution. We treat the axon as a one-dimensional lattice and model the motion of vesicle-bound particles with ordinary differential equations for the expected occupation number at each lattice site. We also assume that each lattice site has a corresponding synapse to which the particle occupying the site can deliver its cargo. In the irreversible case, we use a mean field approximation to recast the original model as a nonlinear partial differential equation reminiscent of the hydrodynamic equations that appear in models of totally asymmetric exclusion processes (TASEP). We find that exclusion effects exacerbate the preferential delivery to proximal synapses when compared to the results of no exclusion obtained in reference [17]. For the reversible delivery case, we allow particles to randomly switch between a motile and stationary state. In contrast to the irreversible case, we also keep track of the motion of particles that are not carrying any vesicles. Hence, the resulting exclusion process has four internal states. The mean field approximation again allows for TASEP-like hydrodynamic equations which, under an adiabatic approximation, can be 45 solved exactly. We find that reversibility in cargo delivery allows for a more homogeneous distribution of vesicles, provided that the presence of a vesicle bound to a motor-cargo complex does not significantly change its speed (hopping rate). TASEP models have been studied extensively by the nonequilibrium statistical physics [6, 8, 22, 28-30, 85, 91, 100, 101] and probability communities [9, 52, 107], and numerous of its applications have been explored. In particular, one can completely describe the time evolution of a system undergoing a TASEP in terms of a chemical master equation and write down differential equations for the expectation of the occupation number for each site in the lattice domain. From these descriptions, mean field approximations and continuum limits can be taken to recast the stochastic dynamics of the TASEP as a nonlinear partial differential equation, which can then be analyzed to obtain information about the steady-state properties [8, 18, 22]. Phase diagrams can then be constructed to describe steady-state behavior in distinct parameter regimes. It is natural to ask: how close to the true dynamics of the TASEP are the dynamics of the approximations? Remarkably, the TASEP is a stochastic process that can be solved exactly using a matrix-product ansatz. Comparison of solutions of the approximate problem to the solutions of the true problem shows good agreement. Hence, it is a convenient framework for analyzing several physical problems. For example, the TASEP is a good representation of traffic flow [4, 31, 128] and the growth of random polymers [2]; partial differential equations have thus been derived that describe these processes. Conversely, we now understand the existence of the chemical master equation underlying several physical problems with local densitydependent velocities [105]. The TASEP is growing in prominence in the theoretical biology community as well [22] and therefore has a natural place in its application to our work. The structure of the chapter is as follows. In section 3.1, we introduce our single-state model of irreversible vesicular transport with exclusion, and show how it maps on to a TASEP. We then turn to the four-state model of reversible vesicular transport with exclusion (section 3.2). Finally, in section 3.3, we briefly relate our model to other models that investigate driven exclusion processes, where internal states are assigned to particles occupying each lattice site, for example [100, 101]. However, it should be noted that in contrast to these other studies, we are not concerned with constructing phase diagrams as a function of model parameters such as the inward and outward fluxes. Rather, we 46 are interested in the particular question of how exclusion effects alter the steady-state distribution of synaptic vesicles. 3.1 Irreversible Vesicular Transport with Exclusion Consider a motor-cargo complex hopping unidirectionally along a one-dimensional track; see Figure 3.1. We represent the track as a lattice of N sites, labeled i = 1, . . . , N, with lattice spacing ε = L/N, where L is the length of the track. For simplicity, we assume that each particle can only carry a single cluster of vesicles, and that we ignore partial delivery of a cluster, that is, it is "all-or-none." In the following, we represent a vesicular cluster by a single vesicle. Each site is either vacant or occupied by a vesicle-bound particle, and the particle can hop to the right if and only if the adjacent site is vacant (hard exclusion). At each site, a particle can irreversibly deliver its vesicle(s) to a synaptic target at a rate K and the corresponding site becomes vacant. In other words, we assume that a motorcargo complex without vesicles does not obstruct the movement of other particles. (This simplification will be removed in our full model; see section 3.2). We specify the state of the site i in terms of the occupation number ni ∈ {0, 1} with ni = 1 if the i-th site is occupied by a vesicle-bound motor-cargo complex and zero otherwise. The hopping rate of a particle is taken to be h. We assume that particles are injected on the left-hand boundary at a rate α, and exit the right-hand boundary at a rate β with 0 < α, β < h. Finally, we assume that each lattice site i 6= 1, N has an associated synaptic target with ci vesicles (taken to be large so that ci is treated as a continuous variable). Within the context of intracellular motor transport, one typically interprets the particle as a single molecular motor and the track as a single microtubular filament, with the fundamental length-scale (lattice spacing) given by a single step of a motor, which is around 10 nm [22]. Here, however, we are interested in the transport of motor-cargo complexes along axons, and the delivery of vesicular cargo to synaptic targets. This means that we are looking at processes occurring on significantly longer length-scales. First, we take a single particle to be a macromolecular complex consisting of multiple motors bound to a cargo. Such a complex could have a size of around 0.1 − 1µm, which is comparable to the size of a synaptic target. Therefore, for concreteness, we take the lattice spacing to be ε = 1µm. Second, the 1D track is now identified with an axon of length L that could 47 synaptic targets h K α β motor-cargo complex vesicle Figure 3.1. Dynamical rules for irreversible vesicular transport: hopping, irreversible exchange of vesicles with synaptic targets, and entry/exit rates. extend for several mm. (For simplicity, we assume that the transfer of motors from one MT to the next along an axon is smooth.) It is important to note that one major simplification of our discrete hopping model is that we are replacing a single continuous run of the motor-cargo complex by a single hop over a distance of ε. We are also assuming that the particle stops at regularly spaced synaptic sites. A more complex, hetereogeneous model would distinguish between the size of the complex, the spacing of synaptic targets, and the fundamental lattice spacing. We are interested in determining macroscopic properties of the above exclusion process, in particular, the steady-state density profiles (average occupancies of each lattice site) and the distribution of synaptic vesicles. The density of motor-cargo complexes is denoted by hni i. Here the angular brackets denote the average with respect to all histories of the stochastic dynamics, which can be interpreted as an ensemble average over a large set of trials starting from the same initial conditions. Away from the boundaries, the dynamics is described by the following system of equations for 1 < i < N: d h ni i = hni−1 (1 − ni )i − hni (1 − ni+1 )i − K hni i. dt (3.1) At the boundaries, we have d h n1 i = −hn1 (1 − n2 )i + αh1 − n1 i, dt dhn N i = hn N −1 (1 − n N )i − βhn N i. dt (3.2) (3.3) Note that we have fixed the unit of time so that the hopping rate h = 1. The number of vesicles at the i-th synaptic target is taken to evolve according to the simple first-order 48 kinetic scheme dci = K hni i − γci , dt (3.4) where γ is a vesicular degradation rate. (As highlighted in section 2.1.1, if we were to neglect degradation of synaptic vesicles, then we would have to impose a maximum capacity of synaptic targets, otherwise ci could become unbounded. This is not an issue for reversible vesicular transport.) Note that if K = 0 (no delivery of vesicles to synaptic targets), then equation (3.1) reduces to the standard totally asymmetric exclusion process (TASEP) [8, 28-30, 66]. On the other hand, if K > 0, then it is equivalent to a limiting case of TASEP with langmuir kinetics [32, 91], in which the motor binding rate is zero We will analyze the above model by using the hydrodynamic approach of Parmeggiani et al. [32, 91]. As is well known, equation (3.1) constitutes a nontrivial many-body problem, since in order to calculate the time evolution of hni i, it is necessary to know the two-point correlations hni−1 ni i. The latter obey dynamical equations involving three-point and fourpoint correlations. Thus, there is an infinite hierarchy of equations of motion. However, progress can be made by using a mean-field approximation and a continuum limit in order to derive a partial differential equation (PDE) for the densities. The mean-field approximation consists of replacing two-point correlations by products of single-site averages: hni n j i = hni ihn j i. (3.5) Next we set x = kε and ρ( x, t) = hnk (t)i. The continuum limit is then defined according to N → ∞ and ε → 0 such that the length of the track L = Nε is fixed. (We fix length scales by setting L = 1.) Taylor expanding ρ( x ± ε, t) in powers of ε, 1 ρ( x ± ε, t) = ρ( x ) ± ε∂ x ρ( x, t) + ε2 ∂ xx ρ( x, t) + O(ε3 ), 2 (3.6) then gives to leading order in ε the following nonlinear PDE: ∂ρ ∂J ( x, t) = −ε − Kρ( x, t), ∂t ∂x (3.7) where J ( x, t) = ρ( x, t)(1 − ρ( x, t)) − ε ∂ρ( x, t) , 2 ∂x (3.8) and the boundary conditions are J (0, t) = α(1 − ρ(0, t)), J (1, t) = βρ(1, t). (3.9) 49 Finally, the continuum limit of equation (3.4) is ∂c( x, t) = Kρ( x, t) − γc( x, t). ∂t 3.1.1 (3.10) Steady-state Analysis We wish to calculate the steady-state distribution of synaptic vesicles, which is given by c( x ) = Kρ( x ) , γ (3.11) with ρ( x ) the solution of the steady-state equation ε Kρ (1 − 2ρ)∂ x ρ − ∂ xx ρ = − . 2 ε (3.12) Following Parmeggiani et al. [91], we drop the O(ε) diffusion term and write the first-order ordinary differential equation (ODE) in the form ∂ x [2ρ( x ) − ln ρ( x )] = K . ε (3.13) The resulting boundary value problem is overdetermined as one still has to satisfy the boundary conditions at x = 0, 1: ρ(0) = α, ρ(1)(1 − ρ(1)) = βρ(1). (3.14) Note that the second boundary condition is satisfied if ρ(1) = 1 − β or ρ(1) = 0. The standard procedure is to separately solve the ODE in the two domains [0, x ) and ( x, 1], imposing the left and right boundary conditions, respectively. The two solutions are matched in an O(ε) neighborhood of some point x0 using a boundary layer. (Within the boundary layer, the density changes rapidly and one can no longer ignore the diffusion term.) This matching also determines the location of x0 . In our particular system, the physically relevant solutions decay (faster than) exponentially from the left-hand boundary x = 0 with some correlation length ξ (see below). Since ξ L, it follows that we can effectively treat the domain as semi-infinite with ρ( x ) → 0 as x → ∞. In particular, the solution is independent of β. 50 Integrating equation (3.13) in the two domains yields the left-end (l) and right-end (r) solutions ρ( x )e−2ρ(x) = Yl,r ( x ), (3.15) with Yl ( x ) = ρ(0)e−Kx/ε−2ρ(0) , Yr ( x ) = ρ(1)e−K(x−1)/ε−2ρ(1) . (3.16) As noted in reference [91], equation (3.15) has an explicit solution expressed in terms of the so-called Lambert W function, 2ρ( x ) = −W (−Y ( x )) with Y ( x ) = 2Yl,r ( x ). The Lambert W function [24] is a multivalued function with two real branches as shown in Figure 3.2. Since ρ( x ) ∈ [0, 1], it follows that  1     − 2 W0 (−Y ( x )), ρ ∈ [0, 0.5], ρ( x ) =     − 1 W−1 (−Y ( x )), ρ ∈ [0.5, 1]. 2 (3.17) In contrast to reference [91], we do not assume that the degradation rate K is O(ε) since this would yield unrealistically slow delivery rates (see below). This means that the left-end function Yl ( x ) decays over a length-scale ξ (in physical units) such that ξ ∼ hL/(KN ). If we take the effective length of the axon to be 10 mm, the lattice spacing to be 1 µm, and the hopping rate to be 0.1 − 1s−1 (based on speeds of motor-cargo complexes [48]), then ξ ∼ K −1 µm with K measured in s−1 . Thus, in order to have correlation lengths comparable to axonal lengths of several mm, we would require delivery rates on the order of k ∼ 10−4 − 10−5 s−1 , whereas measured rates tend to be of the order of a few inverse minutes [45, 74]. Therefore, in contrast to [91], ξ L. Hence, Yl ( x ) ≈ 0 when ξ x < L. Similarly, the right-end function Yr ( x ) grows exponentially over a distance ξ from x = 1. It is clear that the only physically relevant solution when α < 1/2 is ρ( x ) = −W0 (−2Yl ( x ))/2 with ρ(0) = α and ρ(1) = 0. (Since W0 (−Y ) is a monotonically decreasing function of |Y | with W0 (−Y ) → 0 as Y → 0, it follows that the density ρ( x ) also decays over the length-scale within the bulk of the domain.) If α > 1/2, then the left-end solution ρ( x ) = −W1 (−2Yl ( x ))/2 cannot match the right-hand boundary condition, since W−1 (−2Yl ) → ∞ as Yl → 0. Hence, there exists a boundary layer on the left-hand side that matches bl ( x ))/2. Here Y bl ( x ) = ρ(0) = α > 1/2 with a bulk solution of the form ρ( x ) = −W0 (−2Y Ae−Kx/ε , with the constant A determined by matching the solutions in the boundary layer. 51 1 W0(Y) -1/e -0.5 0 0.5 1 1.5 2 Y -1 branch point -2 W-1(Y) -3 -4 -5 Figure 3.2. The real branches W0,−1 (Y ) of the Lambert W function. The main conclusion of the above analysis is that when the delivery of vesicles to synaptic targets is irreversible, with motor-cargo complexes injected at the left-hand side, there is an exponential-like decrease in the distribution of synaptic vesicles along the axon as previously observed in a model without exclusion [17], except that the decay is faster with exclusion. This indicates that exclusion effects exacerbate the preferential delivery of cargo to proximal synapses; see Figure 3.3. A heuristic explanation is that particles move more slowly as they are blocked by exclusion, and will thus be closer to the entrance when they deliver their vesicle. 3.2 Reversible Vesicular Transport with Exclusion We now turn to our full model that combines reversible cargo delivery, exclusion effects and different motile states. As with the simpler advection-diffusion model given by equations (2.5), we now have to keep track of motors with and without vesicular cargo. As with the previous exclusion model (section 3.1), we assume that each particle can only carry a single cluster of vesicles, and that exchange of vesicles is "all-or-none." We also assume that each particle can switch between two states, a motile state (+) and a stationary state (0). When in the stationary state, the particle can reversibly exchange a vesicle with a synaptic target. Again we represent the 1D track as a lattice of N sites, labeled i = 1, . . . , N, with lattice spacing ε = L/N, where L is the length of the track. Each site is either vacant 52 Lambert Exponential 0.4 0.35 0.3 0.25 ρ 0.2 0.15 0.1 0.05 0 0 0.2 0.4 0.6 0.8 1 x [100 μm] Figure 3.3. Comparison of the steady-state solution to equation (3.12) and the decaying exponential seen in reference [17]. Parameter values are β = 0.9, ε = 0.01, and α = 0.4. or occupied by a particle in the motile or stationary state and with or without a vesicle. A motile particle can hop to the right if and only if the adjacent site is vacant (free of any particles). In order to keep track of whether or not a vesicle is bound to a particle, we specify the state of the site i in terms of the occupation numbers ni+,0 ∈ {0, 1} and mi+,0 ∈ {0, 1}. Here ni+,0 = 1 if the i-th site is occupied by a particle in state (+, 0) that is carrying a vesicle, whereas mi+,0 = 1 is the corresponding case when the particle is without a vesicle. The vacancy occupation number χi is then determined by the conservation law χi + ni+ + n0i + mi+ + m0i = 1. (3.18) The hopping rate of a particle is taken to be h if it is carrying a vesicle and by h if it is not. It remains to specify the transition rates between the different internal particle states. First, a particle can switch between the motile and stationary states with rates κ± so that κ− −* (ni+ = 1, n0i = 0) ) − (ni+ = 0, n0i = 1), (3.19) κ+ and κ− −* (mi+ = 1, m0i = 0) ) − (mi+ = 0, m0i = 1). κ+ (3.20) 53 For simplicity, we take the transition rates to be the same whether or not a vesicle is bound to the particle. Second, a vesicle can be reversibly exchanged with a synaptic target according to the rates K± so that K− ci − * (n0i = 0, m0i = 1) − ) − − (n0i = 1, m0i = 0). K+ (3.21) We assume that the number of vesicles ci at the i-th synaptic target is sufficiently large so that it is never depleted. Finally, particles with (without) a bound vesicle are injected on the left-hand boundary at a rate αn (αm ), and exit the right-hand boundary at a rate β. The various processes are illustrated in Figure 3.4. Following along analogous lines to section 3.1, we represent the average with respect to all histories of the stochastic dynamics by angular brackets, and denote the density of particles with (without) a bound vesicle and in the motile state (+) or stationary state (0) by hni+,0 ( a)i (hmi+,0 ( a)i). Away from the boundaries, the dynamics is described by the following system of equations: dhni+ i = hhni+−1 (1 − ni+ − n0i − mi+ − m0i )i dt dhn0i i dt − hhni+ (1 − ni++1 − n0i+1 − mi++1 − m0i+1 )i + κ+ hn0i i − κ− hni+ i, = −κ+ hn0i i + κ− hni+ i + K+ ci hm0i i − K− hn0i i, (3.22) (3.23) and dhmi+ i = hhmi+−1 (1 − ni+ − n0i − mi+ − m0i )i dt dhm0i i dt − hhmi+ (1 − ni++1 − n0i+1 − mi++1 − m0i+1 )i + κ+ hm0i i − κ− hmi+ i, = −κ+ hm0i i + κ− hmi+ i − K+ ci hm0i i + K− hn0i i. (3.24) (3.25) At the boundaries, equations (3.22) and (3.24) become dhn1+ i dt dhn+ Ni dt dhm1+ i dt dhm+ Ni dt = −hhn1+ (1 − n2+ − n02 − m2+ − m02 )i + αn h1 − n1+ − n01 − m1+ − m01 i, (3.26) + 0 + 0 + = hhn+ N −1 (1 − n N − n N − m N − m N )i − β h n N i, (3.27) = −hhm1+ (1 − n2+ − n02 − m2+ − m02 )i + αm h1 − n1+ − n01 − m1+ − m01 i, (3.28) + 0 + 0 + = hhm+ N −1 (1 − n N − n N − m N − m N )i − β h m N i. (3.29) 54 synaptic targets h αn K+ Kβ αm vesicle motile state κ+ κ- stationary state Figure 3.4. Dynamical rules for reversible vesicular transport: hopping, switching between motile and stationary particle states, reversible exchange of vesicles with synaptic targets, and entry/exit rates. Finally, given these densities, the number of vesicles at the i-th synaptic target is taken to evolve according to the simple first-order kinetic scheme dci = K− hn0i i − K+ ci hm0i i. dt 3.2.1 (3.30) Mean-field and Continuum Limit Equations (3.22)-(3.29) constitute a nontrivial many-body problem, since in order to calculate the time evolution of hni+ i, it is necessary to know the two-point correlations hni+−1 ψi i, where ψi ∈ {ni+,0 , mi+,0 } and similarly for hmi+ i. The latter obey dynamical equations involving three-point and four-point correlations. Thus, there is an infinite hierarchy of equations of motion. However, progress can be made by using a mean-field approximation and a continuum limit in order to derive a partial differential equation (PDE) for the densities [32, 91]. The mean-field approximation consists of replacing twopoint correlations by products of single-site averages: hni+ ψj i = hni+ ihψj i, hmi+ ψj i = hmi+ ihψj i. (3.31) ,0 +,0 Next we set x = kε, ρ+,0 ( x, t) = hn+ k ( t )i and σ+,0 ( x, t ) = h mk ( t )i. The continuum limit is then defined according to N → ∞ and ε → 0 such that the length of the track L = Nε is fixed. (We fix length scales by setting L = 1). Taylor expanding ρ+,0 ( x ± ε, t) and σ+,0 ( x ± ε, t) in powers of ε, 1 ρ0 ( x ± ε, t) = ρ0 ( x ) ± ε∂ x ρ0 ( x, t) + ε2 ∂ xx ρ0 ( x, t) + O(ε3 ), 2 (3.32) 55 etc., then gives to leading order in ε the following system of PDEs: ∂Jρ ( x, t) ∂ρ+ = −ε + + κ + ρ0 − κ − ρ + , ∂t ∂x ∂ρ0 = −κ+ ρ0 + κ− ρ+ + K+ cσ0 − K− ρ0 , ∂t (3.33) (3.34) and ∂Jσ ( x, t) ∂σ+ = −ε + + κ+ σ0 − κ− σ+ , ∂t ∂x ∂σ0 = −κ+ σ0 + κ− σ+ − K+ cσ0 + K− ρ0 . ∂t (3.35) (3.36) The currents are Jρ+ = hLρ+ , Jσ+ = hLσ+ , (3.37) where for any function F, LF = (1 − ρ − σ) F − ε [(1 − ρ − σ)∂ x F − F∂ x (1 − ρ − σ)] , 2 (3.38) for ρ = ρ0 + ρ+ and σ = σ0 + σ+ . From equations (3.26)-(3.29), we have the corresponding boundary conditions Jρ+ (0, t) = αn (1 − ρ(0, t) − σ(0, t)), Jσ+ (0, t) = αm (1 − ρ(0, t) − σ (0, t)), (3.39) and Jρ+ (1, t) = βρ+ (1, t), Jσ+ (1, t) = βσ+ (1, t). (3.40) Finally, the continuum limit of equation (3.30) is ∂c( x, t) = K− ρ0 ( x, t) − K+ c( x, t)σ0 ( x, t). ∂t 3.2.2 (3.41) Fast Switching Limit We now make the additional simplification that the rates κ± of switching between the stationary and motile states are much faster than the hopping rate and K± . This is made 56 explicit by performing the rescalings κ± → κ± /δ, where δ is a second small parameter. We can then carry out an adiabatic reduction of equations (3.33) and (3.34) by setting ρ+ ( x, t) = κ+ ρ( x, t) + δw+ ( x, t), κ ρ0 ( x, t) = κ− ρ( x, t) + δw0 ( x, t), κ (3.42) σ+ ( x, t) = κ+ σ( x, t) + δw+ ( x, t), κ σ0 ( x, t) = κ− σ( x, t) + δw0 ( x, t), κ (3.43) and with κ = κ+ + κ− , w0 + w+ = 0, and w0 + w+ = 0. Substituting these expansions into equations (3.33)- (3.36) gives ∂Jw+ ( x, t) κ+ ∂ρ ∂w+ κ+ ∂Jρ ( x, t) +δ = −ε − εδ + κ + w0 − κ − w + , κ ∂t ∂t κ ∂x ∂x ∂w0 κ− ∂ρ +δ = − κ + w0 + κ − w + κ ∂t ∂t κ− + (K+ cσ − K− ρ) + δ (K+ cw0 − K− w0 ) , κ (3.44) (3.45) and ∂Jw+ ( x, t) κ+ ∂σ ∂w+ κ+ ∂Jσ ( x, t) +δ = −ε − εδ + κ + w0 − κ − w + , κ ∂t ∂t κ ∂x ∂x κ− ∂σ ∂w0 +δ = −κ + w0 + κ − w + κ ∂t ∂t κ− − (K+ cσ − K− ρ) − δ (K+ cw0 − K− w0 ) . κ (3.46) (3.47) Here Jρ = hLρ, Jσ = hLσ etc. Adding equations (3.44) and (3.45) yields ∂Jw+ ( x, t) κ− ∂ρ κ+ ∂Jρ ( x, t) = −ε − εδ + (K+ cσ − K− ρ) + δ (K+ cw0 − K− w0 ) , ∂t κ ∂x ∂x κ (3.48) whereas adding equations (3.46) and (3.47) gives, on dropping O(εδ) terms, ∂Jw+ ( x, t) κ− ∂σ κ+ ∂Jσ ( x, t) = −ε − εδ − (K+ cσ − K− ρ) − δ (K+ cw0 − K− w0 ) . ∂t κ ∂x ∂x κ (3.49) Next we substitute for ∂ρ/∂t in equation (3.44) using equation (3.48), substitute for ∂σ/∂t in equation (3.46) using equation (3.49), and introduce the double asymptotic expansions w0 = w0,0 + ∑ i,j,i + j>0 δi ε j w0,ij , w0 = w0,0 + ∑ i,j,i + j>0 δi ε j w0,ij , (3.50) 57 with w+ = −w0 , w+ = −w0 . The lowest order coefficients are w0,0 = κ− κ+ (K+ cσ − K− ρ) , κ2 κ (3.51) and w0,0 = − κ− κ+ (K+ cσ − K− ρ) . κ2 κ (3.52) Hence, equations (3.48) and (3.49) have the leading order form and ∂ρ κ+ ∂Jρ ( x, t) b b− ρ, = −ε + K+ cσ − K ∂t κ ∂x with κ+ ∂Jσ ( x, t) b ∂σ b− ρ, = −ε − K+ cσ + K ∂t κ ∂x (3.53) (3.54) i h b− = κ− K− 1 + δ(K− + cK+ ) κ− κ+ , K κ κ2 i h κ κ κ+ − − b+ = K K+ 1 + δ(K− + cK+ ) 2 . κ κ (3.55) (3.56) Finally, equation (3.41) becomes ∂c( x, t) b− ρ( x, t) − K b+ c( x, t)σ( x, t). =K ∂t (3.57) We note that if h = h, then adding equations (3.53) and (3.54) yields a hydrodynamic equation for the total density of particles φ( x, t) = ρ( x, t) + σ( x, t) identical in form to the totally asymmetric exclusion process (after rescaling): ∂φ( x, τ ) ∂J ( x, τ ) = −ε , ∂τ ∂x (3.58) with J ( x, τ ) = Jρ ( x, τ ) + Jσ ( x, τ ) = φ( x, τ )(1 − φ( x, τ )) − ε ∂φ( x, τ ) , 2 ∂x (3.59) and boundary conditions J (0, t) = α(1 − φ(0, t)), J (1, t) = βφ(1, t), (3.60) The rescalings are τ=ε κ+ t, κ α= κ ( α m + α n ), κ+ h = h = 1. (3.61) 58 3.2.3 Steady-state Analysis We now establish that a uniform, steady-state distribution of synaptic vesicles occurs when h = h = 1. The steady-state equations are c( x ) = b− ρ( x ) K , b+ σ( x ) K (3.62) ε [(1 − φ)∂ x ρ + ρ∂ x φ] = Jρ , 2 ε (1 − φ)σ − [(1 − φ)∂ x σ + σ∂ x φ] = Jσ . 2 (1 − φ ) ρ − (3.63) (3.64) Here Jρ and Jσ are constant nonequilibrium currents for the ρ and σ particles. Adding equations (3.2.3b,c) yields the steady-state version of the TASEP equation (3.58): φ (1 − φ ) − ε dφ = J, 2 dx (3.65) with J = Jρ + Jσ . From the boundary conditions (3.39) and (3.40), it follows that Jρ = καn J, κ+ α Jσ = καm J, κ+ α (3.66) and, hence, equations (3.2.3b,c) have the solution ρ( x ) = καn φ ( x ), κ+ α σ( x) = καm φ ( x ). κ+ α (3.67) Finally, substituting this solution into equation (3.2.3a) yields the constant vesicular distribution c ( x ) = c0 = b− αn K . b+ αm K (3.68) Since both densities σ( x ) and ρ( x ) are proportional to the steady-state solution of the standard TASEP, it is worthwhile briefly recapping the well-known properties of the latter [8, 66]. This will be useful when comparing the corresponding profiles when h 6= h. Setting q = φ − 1/2, the steady-state current equation (3.65) takes the form (after absorbing the factor of 2 into ε) ε dq = v2 − q2 , dx It follows that for v2 > 0 ε Z 1 − J0 . 4 (3.69) dq = x − x0 , (v − q)(v + q) (3.70) v2 = where x0 is an integration constant. Using partial fractions, we find that v+q = e2v(x−x0 )/ε , v−q (3.71) 59 which on rearranging yields the density profile φ( x ) = 1 + v tanh(v( x − x0 )/ε), 2 (3.72) with v ≥ 0. On the other hand, if v2 < 0, then we have ε Z dq | v2 | + q2 = x − x0 . (3.73) Under the change of variables q = cotan(u), we can evaluate the integral and find that φ( x ) = 0.5 + |v|cotan(|v|( x − x0 )/ε). (3.74) The two unknown parameters J0 , x0 can be determined in terms of α, β by imposing the boundary conditions at x = 0, L. As is well known, three distinct phases can be identified [8, 66] (see Figure 3.5(d)): 1. A low-density phase in which the bulk density is smaller than 1/2, x0 ≈ 1 and v2 > 0. Since ε 1, we see from equation (3.72) that φ( x ) ≈ 0.5 − v for all x < x0 . In particular, at the left-hand boundary α(0.5 + v) = J0 , which can be rewritten as v = J0 /α − 0.5. Squaring both sides and using the definition of v gives, to lowest order in ε, φ(0) = α, J0 = α(1 − α), α < 1/2. (3.75) The other boundary condition becomes β= J0 J0 > = α. 0.5 + v tanh(v( L − x0 )/ε) 0.5 + v (3.76) In order to satisfy this boundary condition, there is an ε-wide boundary layer at x = L with L − x0 = O(ε). 2. A high-density phase in which the bulk density is larger than 1/2 and x0 ≈ 0. Hence, φ( x ) ≈ 0.5 + v in the bulk of the domain and at the right-hand boundary we have β(0.5 + v) = J0 . Following along similar lines to the low-density case, we deduce that φ( L) = 1 − β, J0 = β(1 − β), β < 1/2, (3.77) and β < α. There is now a boundary layer around x = 0 in order to match the rate α. The two phases coexist along the line α = β < 1/2. 60 (a) φ 1 (b) 1 0.9 0.9 0.8 0.8 0.7 0.7 0.6 φ 0.5 0.6 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0.2 0.4 0.6 0.8 1 0 0 0.2 0.4 x [100 μm] (c) 0.6 0.8 1 x [100 μm] 1 (d) 0.9 1 J0 =1/4 J0 =α(1−α) 0.8 0.7 LD 0.6 φ 0.5 MC β 0.5 0.4 J0=β(1−β) 0.3 0.2 HD 0.1 0 0 0.2 0.4 0.6 0.8 1 0 x [100 μm] 0.5 α 1 Figure 3.5. Steady-state solutions of total density φ in the different phases with ε = 0.01 and L = 1. (a) Plot of φ in the HD phase for α = 0.9 and β = 0.3. (b) Plot of φ in the LD phase for α = 0.1 and β = 0.6. (c) Plot of φ in the MC phase for α = β = 0.7. (d) Mean-field phase diagram for the TASEP showing the regions of α, β parameter space where the low-density (LD), high-density (HD) and maximal-current (MC) phases exist. 3. A maximal current phase. In the region α > 1/2, β > 1/2, we require J0 > 1/4 so that v2 < 0. It turns out that the current takes the form J0 = 0.25 + O(ε2 /L2 ), that is, it is very close to the maximal value of function φ(1 − φ). This follows from the observation that the solution (3.74) will blow up unless 0 < |v|( x − x0 )/ε < π for all x ∈ [0, L]. This implies that x0 = −O(ε) and |v| < πε/L. Under these conditions, equation (3.74) ensures that φ( x ) ≈ 0.5 in the bulk of the domain. The precise values of v and x0 are then adjusted so that the boundary conditions at x = 0, L are satisfied: φ(0) = 1 − 1/(4α) > 0.5 and φ( L) = 1/(4β) < 0.5. Also note that away from the left-hand boundary, we have cotan(|v|( x − x0 )/ε) ≈ ε/(|v| x ) so that φ( x ) ∼ 0.5 + ε/x. (3.78) In deriving equations (3.53) and (3.54), we first adopted the mean-field approximation used to study TASEP models with single internal states [32, 91], and then carried out an 61 1 a.) 0.9 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 1 0.2 0.4 0.6 0.8 0 1 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 c 0.9 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0.2 0.4 0.6 0.8 1 1 c.) 0.9 0 b.) 0.9 σ ρ 0.8 φ 1 κ = 10 + κ+ = 1 κ + = .1 0.2 0.4 0.6 x [100 µ m] 0.8 1 0 d.) 0 0.2 0.4 0.6 0.8 1 x [100 µm] Figure 3.6. Effect of slowing down the switching rates between motile and immotile states on concentration profiles when TASEP limit is in a maximum current phase. Plots of (a) ρ, (b) σ (c) total motor density φ, and (d) and synaptic vesicle density c for various switching rates κ− = κ+ . Other parameter values are αn = αm = 0.8, β = 0.8, K± = 0.5, h = h = 1, and N = 100. adiabatic approximation in the fast switching limit. If these approximations are valid, then we expect numerical simulations of the full stochastic model to generate a total motor density profile φ that converges to the classical TASEP density in the limit κ± → ∞ for h = h. This is indeed found to be the case as illustrated in Figure 3.6(c). We can see that the profile of φ for fast switching in the maximal current parameter regime resembles the profile for the classic TASEP model. However, as the switching slows down, the profile deviates from the TASEP curve. Nevertheless, this does not have a significant effect on the distribution of synaptic vesicles, since c is still approximately uniform; see Figure 3.6. Interestingly, it has been shown in references [92] and [93] that standard mean field theory can break down for a model in which particles switch between motile and stationary states, due to statistical correlations between motile and stationary occupation 62 numbers. Numerically , we find that this does not present a problem for our particular model when the system operates in a regime where the switching rates κ± between the motile and immotile states of the motors are fast compared to the hopping rate h and rates of exchange of vesicles between motors and synapses K± . Note that Figure 3.6 and subsequent numerically generated figures are generated using a continuous-time Monte Carlo algorithm based on the Gillespie algorithm [35] and the dynamical rules elucidated in Figure 3.4. Individual particles carrying cargo that are bound to a microtubule can move to the adjacent site at a rate h provided the adjacent site is unoccupied. Particles not carrying cargo but bound to a microtubule can move to the adjacent site provided it is empty at a rate h. Individual particles may bind and unbind from a microtubule at the rates κ± and particles unbound from microtubules may deliver vesicles at a rate K− or recover them at a rate K+ . We collect statistics from the system once it has reached steady state. To ensure it has reached steady state, we neglect the first 108 steps and collect statistics on the subsequent 108 steps. 3.2.4 Disparity in Hopping Rates In the case where h 6= h, adding together equations (3.53) and (3.54) yields h ∂J ∂φ ∂Jσ i ρ , = −ε + ∂t ∂x ∂x (3.79) which cannot be easily analyzed. Nevertheless, the time evolution of the system can be understood by performing Monte Carlo simulations of the full stochastic model as summarized above. We find that the value of H ≡ h − h alters the nature of the distribution of vesicles along the axon. This is illustrated in Figure 3.7, Figure 3.8, and Figure 3.9, which correspond respectively to the LD, HD, and MC phases for φ in the limit h = h = 1. In each figure, we plot the density profiles of ρ, σ, φ, and c for various hopping rates h < h = 1. It can be seen that in each case, as h decreases (H increases), the distribution c of synaptic vesicles along the axon develops an exponential-like decay with respect to x. This reflects the fact that the ratio ρ( x )/σ ( x ) is no longer x-independent. When h = h, the synaptic vesicle concentration is uniform, c( x ) = 1. We conclude that achieving synaptic democracy is also dependent on the motility of the motor-cargo complexes relative to the motility of the particles without vesicles. In all the stochastic simulations, we take h, the hopping rate of vesicle-bound particles, to be at most h, the hopping rate of particles without vesicles, 63 1 1 a.) 0.9 TASEP h = 0.75 h = 0.5 h = .01 0.8 0.7 0.7 0.6 0.6 0.5 0.5 σ ρ 0.8 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0.2 0.4 0.6 0.8 0 1 1 0.8 0.7 0.7 0.6 0.6 0.5 0.5 c φ 0.2 0.4 0.6 0.8 1 d.) 0.9 0.8 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 1 c.) 0.9 b.) 0.9 0 0.2 0.4 0.6 x [100 µ m] 0.8 1 0 0 0.2 0.4 0.6 0.8 1 x [100 µ m] Figure 3.7. Effect of disparity in hopping rates on concentration profiles when TASEP limit is in a low-density phase. Plots of (a) ρ, (b) σ (c) total motor density φ, and (d) and synaptic vesicle density c for various hopping rates h ≤ h = 1. Other parameter values are αm = αn = 0.4, β = 0.7, K± = 0.5, κ± = 10, and N = 100. which corresponds to the intuition that the former would naturally move slower than the latter due to the added load. Hence, there is a correlation between the value of h and the specific type of cargo being delivered. If, for example, the cargo of a motor is too large, then we expect h h, and the distribution of the given cargo along the axon may not be uniform. On the other hand, if the cargo is relatively small, then h ≈ h and synaptic democracy can be achieved. Analogous results were found in reference [17] for the simpler model without exclusion. 3.3 Relationship with Other Exclusion Process Models Equations (3.53) and (3.54) closely resemble the hydrodynamic equations that arise in modeling processes that account for exclusion effects as well as internal states. For example, Reichenbach et al. [100, 101] allow for particles in each lattice site to exist in one 64 1 1 a.) TASEP h = 0.75 h = 0.5 h = .01 0.9 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 1 0.2 0.4 0.6 0.8 0 1 0.2 0.4 0.6 0.8 1 d.) 0.9 0.8 0.7 0.7 0.6 0.6 0.5 0.5 c 0.8 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 1 c.) 0.9 φ 0.8 0.7 σ ρ 0.8 b.) 0.9 0 0.2 0.4 0.6 0.8 1 0 x [100 µ m] 0 0.2 0.4 0.6 0.8 1 x [100 µ m] Figure 3.8. Effect of disparity in hopping rates on concentration profiles when TASEP limit is in a high-density phase. Plots of (a) ρ, (b) σ (c) total motor density φ, and (d) and synaptic vesicle density c for various hopping rates h ≤ h = 1. Other parameter values are αm = αn = 0.9, β = 0.1, K± = 0.5,κ± = 10, and N = 100. of two internal "spin" states; see Figure 3.10. Particles with opposite spins can occupy the same lattice point and can move to the next lattice site at a prescribed rate provided the adjacent site is not already occupied by another particle of the same spin state. Hence, each particle respects the Pauli exclusion principle. Another common interpretation for these internal states is that of a car traveling on one lane of a two-lane highway. In this context, each lattice site corresponds to a segment of the highway, and thus can be occupied by two cars so long as they are not on the same lane. In either of the interpretations, particles are allowed to switch states provided they are alone in occupying a given site. Note that the effects of exclusion on collective vesicle transport has also been analyzed by Muhuri and Pagonabarraga. They consider the case of bidirectional transport in which particles can reverse direction and reversibly bind to the filament [86]. However, the authors do not separately model vesicles and molecular motors. 65 1 1 a.) 0.9 TASEP h = 0.75 h = 0.5 h = .01 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0.2 0.4 0.6 0.8 0 1 c.) 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 c 0.9 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0.2 0.4 0.6 0.8 1 1 1 φ 0.8 0.7 σ ρ 0.8 b.) 0.9 0 0.2 0.4 0.6 x [100 µ m] 0.8 1 0 d.) 0 0.2 0.4 0.6 0.8 1 x [100 µm] Figure 3.9. Effect of disparity in hopping rates on concentration profiles when TASEP limit is in a maximum current phase. Plots of (a) ρ, (b) σ (c) total motor density φ, and (d) and synaptic vesicle density c for various hopping rates h ≤ h = 1. Other parameter values are αm = αn = 0.9, β = 0.7, K± = 0.5, κ± = 10, and N = 100. In our work, we provide a new biophysical example of internal states within the context of exclusion processes. The full model without the application of the adiabatic approximation consists of particles in one of four internal states: (i) a motile particle bound to the track and carrying a vesicle, (ii) a motile particle bound to the track without a vesicle, (iii) a stationary particle unbound from the track but carrying a vesicle, and (iv) a stationary particle unbound from the track without a vesicle. One important difference between the spin and traffic models and ours lies in the definitions of the currents in each model. In spite of the existence of two internal states in the two-lane traffic traffic and spin models, the currents are nevertheless the same as seen in standard TASEP models. That is, if ρi ( x, t) is the density of a particle in the i-th internal state, its current is given by, for example, an expression of the form ρi (1 − ρi ). This arises from the fact that double occupation of a single lattice site is allowed provided each particle exists in a different internal state. 66 h h α κ β β α Figure 3.10. Dynamical rules for an exclusion model with two internal spin states [100, 101]. Particles in up (down) states enter with rates α↑ (α↓ ), move unidirectionally to the right with hopping rate h, flip spin state at a rate κ, and leave the system at rates β↑ (β↓ ). Pauli's exclusion principle holds at every lattice site. In our model, currents take a more restrictive form, since a motor can only hop to the adjacent site if it is completely unoccupied. Hence, the currents in our model have the form shown in equation (3.38). Differences in the currents persist when we use an adiabatic approximation to reduce the full model to a model with two internal states (particles with or without a vesicle). 3.4 Discussion In this chapter, we investigated the biophysical machinery involved in maintaining synaptic democracy in axons. In particular, we generalized the results found in reference [17] by examining the effects of exclusion on the distribution of synaptic vesicles along an axon. For both the irreversible and reversible delivery cases, we modeled the dynamics of motor-cargo complexes in terms of the equations of motion for the average occupation numbers at each site on a 1D lattice. By invoking the mean field approximation, we derived a system of hydrodynamic equations which were used to determine the steady-state distributions of both motor-cargo complexes and synaptic vesicles. In the irreversible case, we found that exclusion exacerbates the preferential delivery of vesicles to synaptic sites near the soma. In the reversible case, we performed an adiabatic approximation on the system of hydrodynamic equations by assuming that switching between internal states is fast compared to ballistic dynamics. We found that the steady-state distribution of vesicles is now approximately uniform, provided that the speed of a particle is only weakly dependent on whether or not it is carrying a vesicle. There are a number of possible extensions of our work on exclusion processes and 67 active transport. One assumption we made was that each motor-cargo complex only carried one SVP. It would be interesting to see what happens if we relax this condition and allow each complex to carry a finite number of SVPs. Some work has been done in this area using aggregation theory and a modified version of the well-known Becker-Doring equations [13], but exclusion effects were not taken into account. Another important generalization would be to investigate what happens when we allow for bidirectional motor transport (see discussion in section 2.6). Yet another extension of our work would be to explore the effects of heterogeneity where, for example, the distribution and size of synaptic targets are not uniform. Finally, we hope to investigate the impact of exclusion effects on other biological processes that involve axonal transport. Examples include the models of motor-based length control presented in Chapters 4 and 5, which do not take into account exclusion effects. CHAPTER 4 FLAGELLAR LENGTH CONTROL A fundamental question in cell biology is how the sizes of subcellular stuctures are determined in order to scale with the size of the cell and with physiological requirements. It appears that self-organizing processes together with physical constraints play a major role in controlling organelle size [97]. At least three distinct control mechanisms have been identified [83]. I. Molecular rulers. In the case of linear structures such as filaments, size control can be achieved by a molecular ruler protein, whose length is equal to the desired length of the growing structure. One classical example is the length of the λ-phage tail, which is determined by the size of the gene H product (gpH) [58]. During assembly of the tail, gpH is attached to the growing end in a folded state, and protects the growing end from the terminator gene product U (gpU). As the tail elongates, gpH stretches such that when it is fully extended, further growth exposes the tail to the action of gpU; see Figure 4.1. II. Quantal synthesis. Size could be controlled by synthesizing exactly enough material to build a structure of the appropriate size - a process known as quantal synthesis. For example, precursor protein levels are known to affect the length of flagella in the unicellular green alga Chlamydomonas reinhardtii [73], and the length of sea urchin cilia is correlated with the concentration of the protein tektin [113]. One prediction of the quantal synthesis model is that doubling the number of flagella should halve their length. However, studies of Chlamydomonas mutants indicate a much weaker dependence of length on the number of flagella, suggesting that there is an additional length-controlling mechanism involving dynamic balance [81]; see below. III. Dynamic balance. Dynamic structures are constantly turning over so that in order for them to maintain a fixed size, there must be a balance between the rates of assembly and disassembly. If these rates depend on the size in an appropriate way, then there 69 gpU gpH Figure 4.1. Schematic of a molecular ruler. Redrawn from [12]. will be a unique balance point that stabilizes the size of the organelle. For example, eukaryotic flagellar microtubules undergo continuous assembly and disassembly at their tips, in which a constant rate of disassembly is balanced by a length-dependent rate of assembly due to a fixed number of molecular motors transporting from the cell body, leading to a fixed flagellar length [81, 82, 108]. An analogous dynamic balance mechanism is thought to control the length of actin-based structures, such as the stereocilia of the inner ear [96, 103]. Here actin filaments constantly treadmill back towards the cell body, with disassembly at the base balanced by assembly at the tip. The latter depends on the diffusion of actin monomers to the tip, which results in a length-dependent rate of assembly. It has also been suggested that a diffusion-induced length dependence of the assembly rate plays a role in the control of the hook length in bacterial flagella [59]. A different balance mechanism appears to control the length of microtubules in yeast, where kinesin motors move processively to the microtubule tips where they catalyze disassembly. Longer microtubules recruit more kinesin motors from the cytoplasm, which results in a length-dependent rate of disassembly. When this is combined with a length-independent rate of assembly, a unique steady-state microtubule length is obtained [39, 47, 53, 99, 100, 122]. A related mechanism involves the modulation of microtubular dynamic instabilities, that is, the catastrophe frequency [50, 114]. One class of organism where such a dynamic mechanism may occur is eukaryotic flagella [76, 77, 81, 82, 108, 126]. These are microtubule-based structures that extend to about 10 µm from the cell and are surrounded by an extension of the plasma membrane. They are at least an order of magnitude longer than bacterial flagella. Flagellar length control is a particularly convenient system for studying organelle size regulation, since a flagellum 70 can be treated as a one-dimensional structure whose size is characterized by a single length variable. The length of a eukaryotic flagellum is important for proper cell motility, and a number of human diseases appear to be correlated with abnormal length flagella. Radioactive pulse-labeling has been used to measure protein turnover in eukaryotic flagella. Such measurements have established that turnover of tubulin occurs at the + end of flagellar microtubules, and that the assembly (rather than disassembly) part of the turnover is mediated by intraflagellar transport (IFT). This is a motor-assisted motility within flagella in which large protein complexes move from one end of the flagellum to the other [108]. Particles of various size travel to the flagellar tip (anterograde transport) at around 1.5 µm/s, and smaller particles return from the tip (retrograde transport) at around 2.5 µm/s after dropping off their cargo of assembly proteins at the + end. A schematic diagram of IFT transport is shown in Figure 4.2. One suggested mechanism for flagellar length control is based on the idea that a fixed number of IFT particles is present inside the flagellum [81, 82]. As the flagellum grows in length, each IFT particle has to travel a longer distance to deliver tubulin at the tip of the flagellum, resulting in a balance between assembly and disassembly at a critical flagellar length. More specifically, if a fixed number of transport complexes N move at a fixed mean speed v, then the rate of transport and assembly should decrease inversely with the flagellar length L. On the other hand, measurements of the rate of flagellar shrinkage when IFT is blocked indicate that the rate of disassembly is length-independent. This has motivated the following simple deterministic model for length control [82]: L(t) vmicrotubule degradation + V _ IFT particle v+ cargo insertion Figure 4.2. Schematic diagram of intraflagellar transport (IFT), in which IFT particles travel with speed v± to the ± end of a flagellum. When an IFT particle reaches the + end, it releases its cargo of protein precursors that contribute to the assembly of the flagellum. Disassembly occurs independently of IFT transport at a speed V. 71 avN dL = − V, dt 2L (4.1) where a is the size of the precursor protein transported by each IFT particle and V is the speed of disassembly. Equation (4.1) has a unique stable equilibrium given by L∗ = avN/2V. Using the experimentally-based values N = 10, v = 2 µm/s, L∗ = 10 µm and V = 0.01 µm/s, the effective precursor protein size is estimated to be a ≈ 10 nm (a stochastic version of a model for intraflagellar transport has also been developed using the theory of continuous time random walks [11]). However, recent photobleaching studies have shown that the influx of IFT particles into the flagellum is regulated [76]. Trains of IFTs enter the flagellum through the flagellar pore, a membrane-spanning structure at the base of the flagellum that may be homologous to a nuclear pore. (The latter regulates the exchange of macromolecules between a cell's cytoplasm and the nucleus.) There is also a microtubule-organizing center known as the basal body, which anchors the flagellar microtubules at the plasma membrane and integrates them with the cytoplasmic microtubules. IFT proteins dock around the basal body and assemble into trains prior to entering the flagellum [27]. It appears that the rate at which IFTs enter the flagellum depends on the amount of docked IFTs in the basal body, with faster growing flagella having more localized IFTs [77, 126]. This suggests that there is some length-dependent mechanism for regulating the accumulation of IFT particles at the basal body (and possibly the loading of cargo to docked IFTs [126]). Ludington et al. [76] considered several different mathematical models of IFT regulation, based on the idea that cell signaling within the flagellum results in a length-dependent binding rate of IFTs within the basal body, and compared the models with experimental photobleaching data on the variation of IFT numbers as a function of length. In this chapter, we develop a stochastic model of flagellar length control that incorporates IFT regulation along the lines of Ludington et al. [76]. In particular, we will assume that the binding rate of IFTs is length-dependent and that the rate of flagellar growth is much slower than any other dynamical process (adiabatic approximation). We then take into account two distinct sources of stochasticity. The first is given by fluctuations in the number of bound IFTs within the basal body due to the random nature of binding-unbinding chemical reactions. According to the law of large numbers, these 72 √ fluctuations will vary as 1/ M where M is the total number of binding sites. Fluctuations in the number of bound IFTs will then result in fluctuations in the rate at which IFT particles are injected into the flagellum. However, a second source of noise arises even when such fluctuations are neglected (M → ∞). That is, even for constant injection rates, we expect intrinsic fluctuations in the injection times. Assuming for simplicity, that the injection times are exponentially distributed, the process of IFT injection for constant rates can be modeled as a Poisson process. It immediately follows that if fluctuations in the number of bound IFTs are also taken into account, then the resulting model of IFT length control is given by a doubly stochastic Poisson process (DSPP). DSPPs were first introduced by Cox [25] as a generalization of an inhomogeneous Poisson process, in which the time-dependent transition rate depends on a second, independent stochastic process. The general theory of DSPPs was subsequently developed by Grandell [41]. Example applications include photon and electron detection [104], occurrences of credit events in finance [71], and neural coding [7, 67, 110]. We use the theory of DSPPs to analyze our stochastic model of flagellar length control, and determine how fluctuations in the IFT flux and flagellar number vary with length. The organization of the chapter is as follows. In section 4.1 we introduce our model of flagellar length control and the two major stochastic components, namely, IFT particle binding/unbinding within the basal body and the Poisson-like injection of IFTs into the flagellum. We also present one of the examples of IFT binding regulation considered by Ludington et al. [76]. However, the subsequent analysis is not restricted to any specific regulatory mechanism. In section 4.2, we analyze the statistics of IFTs within the flagellum as a function of flagellar length under the assumption that the number of binding sites is large. The case of a small number of binding sites is considered in section 4.3. 4.1 Stochastic Model of Flagellar Length Control We begin by briefly describing the deterministic model of flagellar length control introduced in reference [76]; see Figure 4.3. Consider a one-dimensional flagellum of length L with the basal body at x = 0 and the tip at x = L. Suppose that there are M binding sites for IFT particles in the basal body, and the concentration of IFTs within the cytoplasm is B. Denote the binding and unbinding rates by k + and k − , respectively. Assuming that M 73 IFT particle basal binding site k+ kflagellum basal body Figure 4.3. Schematic diagram of the basic model. IFT particles (filled circles) can undergo binding/unbinding reactions with M sites (filled rectangles) in the basal body at rates k ± . The number of bound IFTs determines the rate at which IFTs are injected into the flagellum. Once in the flagellum, IFTs are actively transported to the tip, where they deliver their cargo and are then transported back to the basal body along the lines shown in Figure 4.2. Some signaling mechanism within the flagellum (not shown) results in the binding rate k + being dependent on the flagellar length L, resulting in a length-dependent IFT flux regulation. is sufficiently large, the kinetic equation for the number m(t) of bound IFTs at time t is dm = k + B[ M − m(t)] − k − m(t), dt (4.2) which has the steady-state solution m∗ = MX ∗ , X∗ = k+ B . k+ B + k− (4.3) Now suppose that there is some signaling mechanism within the flagellum such that the binding rate is a decreasing function of length L, and set k + = k+ C0 ( L). We will give one example of such a signaling mechanism in section 4.2.3; see also reference [76]. Under the adiabatic approximation that the growth-rate of the flagellum is much slower than the various kinetic processes, we can still treat M∗ as a constant with m∗ = m∗ ( L) ≡ k+ C0 ( L) B M. k+ C0 ( L) B + k − (4.4) The rate of injection of IFTs into the flagellum is then taken to be λ0 = ηm∗ ( L), which means that the influx is a monotonically decreasing function of L. The critical flagellar length is then determined by the balance between the influx and the length-independent 74 rate of disassembly, along analogous lines to equation (4.1). The number of binding sites M ranges from 10 − 1000, whereas fits with experimental data suggest that k+ B/k − ∼ 10µm [76]. 4.1.1 IFT Injection as a Poisson Process The first level of stochasticity occurs if we take the injection times of the IFT particles to be exponentially distributed with rate λ0 = ηm∗ . The number N (t) of particles injected into flagellum over the interval [0, t] is then given by a Poisson process . That is, setting Pn (t) = P( N (t) = n| N (0) = 0), we have dPn = λ0 [ Pn−1 (t) − Pn (t)], dt (4.5) which has the solution Pn (t) = ( λ 0 t ) n − λ0 t e , n! (4.6) It immediately follows that h N (t)i = λ0 t, var[ N (t)] = λ0 t. (4.7) Suppose that each injected particle remains in the flagellum a time T = 2L/v + τ before being removed, where v is the arithmetic mean of the anterograde and retrograde speeds of each IFT particle, and τ is the time spent at the tip. We will take τ = 1s and v = 2µm/s. It follows that the number of particles in the flagellum at time t, t > T is F ( t ) = N ( t ) − N ( t − T ). (4.8) In particular, h F (t)i = λ0 (t − (t − T )) = λ0 T = η k+ C0 ( L) BM (2L/v + τ ). k+ C0 ( L) B + k − (4.9) This yields an expression for how the total IFT protein in the flagellum varies with length L, which can then be compared to photobleaching data along the lines of reference [76]. 75 We can also estimate the variance in the number of IFT particles according to the following: var[ F (t)] = h N (t)2 + N (t − T )2 − 2N (t) N (t − T )i − λ20 T 2 = var[ N (t)] + var[ N (t − T )] − λ20 T 2 + λ20 t2 + λ20 (t − T )2 − 2h N (t) N (t − T )i = λ0 (2t − T ) + 2λ20 t(t − T ) − 2[λ20 (t − T )t + λ0 (t − T )] = λ0 T. (4.10) We have used the fact that for a Poisson process with rate λ0 , the autocorrelation function is h N (t) N (t − T )i = h[ N (t) − N (t − T )] N (t − T )i + h N (t − T )2 i = h[ N (t) − N (t − T )]ih N (t − T )i + var[ N (t − T )] + h N (t − T )i2 = λ20 (t − T )t + λ0 (t − T ). 4.1.2 (4.11) Stochastic Model of IFT Docking at the Basal Body The above analysis ignores fluctuations of m(t) and the fact that when a particle is injected, it leaves a vacant binding site. For the moment, let us ignore the latter effect by assuming that the number of bound sites is much greater than one. (The case of M = O(1) will be considered in section 4.3.) The binding of IFTs is then independent of the Poisson process (but not vice versa). Let Q(m, t) denote the probability that m out of M binding sites are bound by IFTs at time t. The corresponding master equation is dQ(m, t) = k + B( M − m + 1) Q(m − 1, t) + k − (m + 1) Q(m + 1, t) dt − [k + B( M − m) + k − m] Q(m, t), (4.12) with Q(−1, t) = Q( N + 1, t) ≡ 0. This can be rewritten in the form of the birth-death master equation d Q(m, t) = ω+ (m − 1) Q(m − 1, t) + ω− (m + 1) Q(m + 1, t) dt (4.13) − [ω+ (m) + ω− (m)] Q(m, t). with transition rates ω+ (m) = ( M − m)k + B, ω− (m) = mk − . (4.14) 76 A standard calculation yields the steady-state solution Qs (m) of the master equation (4.13) [12, 34]. First, note that it satisfies J (m) = J (m + 1) with J ( m ) = ω − ( m ) Q s ( m ) − ω + ( m − 1) Q s ( m − 1). (4.15) Using the fact that m is a nonnegative integer, that is, Qs (m) = 0 for m < 0, it follows that J (m) = 0 for all m. Hence, by iteration, m ω + ( k − 1) , ω− ( k ) k =1 Q s ( m ) = Q s (0) ∏ with M Q s (0) = m ω + ( k − 1) 1+ ∑ ∏ ω− ( k ) m =1 k =1 ! −1 (4.16) . (4.17) M! . m!( M − m)! (4.18) In the particular case of the transition rates (4.14), we have k+ B Q s ( m ) = Q s (0) k− m After calculating Qs (0), we obtain the binomial distribution Qs (m) = M−m (k + B)m k − M! M! = ( X ∗ ) m (1 − X ∗ ) M − m , M (k + B + k − ) m!( M − m)! m!( M − m)! (4.19) where X ∗ is given by equation (4.3). Using standard formulae for the moments of the binomial distribution, we find that, at steady state, the mean number of bound IFTs at the basal body is hmi = MX ∗ = m∗ , (4.20) where m∗ is the fixed point solution of the kinetic equation (4.4). Similarly, the steady-state variance is with p Var[m] = MX ∗ (1 − X ∗ ), (4.21) √ Var[m]/hmi ∼ 1/ M. Hence, in the large-M limit, we can simply treat the number of bound IFTs as a constant m∗ (for fixed L). The injection of IFTs is then given by a homogeneous Poisson process with rate λ0 = ηm∗ . However, if the total number of binding sites takes intermediate values, M ∼ 100 − 1000 [76], we should really treat m(t) as a stochastic variable evolving according to the birth-death master equation (4.12) and set λ = ηm(t). It follows that the process of IFT injection into the flagellum is described by a so-called doubly stochastic Poisson process (DSPP); see section 4.2.1. 77 In order to facilitate later calculations, we will carry out a system-size expansion of the master equation (4.12) for intermediate values of M [12, 34]. First, introduce the rescaled variable x = m/M and corresponding transition rates MΩ± ( x ) = ω± ( Mx ). Equation (4.12) can then be rewritten in the form dΠ( x, t) = M[Ω+ ( x − 1/M)Π( x − 1/M, t) + Ω− ( x + 1/M)Π( x + 1/M, t) dt − (Ω+ ( x ) + Ω− ( x ))Π( x, t)], (4.22) where Π( x, t) = Q( Mx, t). Treating x as a continuous variable and Taylor expanding terms on the right-hand side to second order in M−1 leads to the Fokker-Planck (FP) equation ∂P( x, t) ∂ 1 = − [ A( x ) P( x, t)] + ∂t ∂x 2M M ∂2 ∑ ∂x2 [ B(x) P(x, t)] , (4.23) k =1 with A( x ) = Ω+ ( x ) − Ω− ( x ) = (1 − x )k + B − k − x, (4.24) B( x ) = Ω+ ( x ) + Ω− ( x ) = (1 − x )k + B + k − x. (4.25) The solution to the FP equation (4.23) determines the probability density function for a corresponding Ito stochastic process X (t), which evolves according to the stochastic differential equation (SDE) dX = A( X )dt + r B( X ) dW (t). M (4.26) Here W (t) denotes an independent Wiener process such that hW (t)i = 0, hW (t)W (s)i = min(t, s). (4.27) We now make the approximation λ(t) = η MX (t). (Certain care must be taken, however, since there is a nonzero probability that X (t) becomes negative. We will assume that this does not cause problems for sufficiently large M.) 4.1.3 RanGTP Model of IFT Flux Regulation The one remaining component of the model is the specification of the length-dependent function C0 ( L) of the IFT binding rate k + . Ludington et al. [76] considered several different signaling mechanisms for generating this length dependence. For the sake of illustration, 78 we will consider one of the models that fit particularly well with photobleaching data. It is a diffusion-based model of RanGTP concentration gradient formation. RanGTP is a small enzyme that is known to play an important role in regulating nuclear transport through the nuclear pore complex, and it is hypothesized that RanGTP plays an analogous role in regulating IFT particle influx. In particular, a decrease in RanGTP concentration at the basal body as cell length increases leads to a reduction in IFT particle influx. Suppose that RanGTP is produced at a rate σ at the tip (x = L), resulting in a concentration gradient; see Figure 4.4. Assume that cytoplasmic RanGTP concentration is negligible and κ is the flow rate through the pore at x = 0. Then the RanGTP concentration per unit volume C ( x, t) evolves as ∂2 C ∂C = D 2 − γC, ∂t ∂x x ∈ [0, L], (4.28) where γ is a degradation rate. The boundary conditions are D ∂C = κC, ∂x x = 0; D ∂C = σ, ∂x x = L. (4.29) Integrating equation (4.28) with respect to x and using the boundary conditions gives dR = σ − κC (0, t) − γR, dt (4.30) where R(t) is the total number of RanGTP molecules per unit area, R(t) = Z L C ( x, t)dx. 0 If we assume that diffusion is fast so that the characteristic length (4.31) p D/γ L, then C ( x, t) is approximately uniform and we can take C (0, t) ≈ R(t)/L. Therefore, dR R = σ − κ − γR. dt L (4.32) Equation (4.32) has the steady-state solution R= σL , γL + κ (4.33) so that the concentration at the basal pore is C0 = C0 ( L) = σ . γL + κ (4.34) Typical values of the parameters are [76] σ ∼ 5 − 20/s, κ ∼ 5 − 25µm/s, γ ∼ 10 − 400/s. Unless stated otherwise, we will take γ/σ = 4 and κ/σ = 1µm. (4.35) 79 k+ κ k- σ RanGTP x=0 flagellum x=L basal body Figure 4.4. Schematic diagram of RanGTP concentration gradient model of IFT flux regulation. A source of RanGTP at the tip of the flagellum sets up a concentration gradient along the flagellum resulting in a length-dependent concentration of RanGTP in the basal body. This in turn regulates the binding rate of IFTs to sites in the basal body. 4.2 Analysis of Model Using the Theory of Doubly Stochastic Poisson Processes In this section, we use the theory of DSPPs to analyze our model of flagellar length control. 4.2.1 Doubly Stochastic Poisson Process Combining the two sources of noise outlined in sections 4.1.1 and 4.1.2, the homogeneous Poisson process given by equation (4.5) becomes a DSPP. That is, { N (t), t ≥ 0} is a counting process with positive intensity λ( X (t)) = ηMX (t), which depends on a second independent stochastic process { X (t), t ≥ 0} with X (t) the fraction of bound IFTs in the basal body. The latter evolves according to the SDE (4.26). For a given realization of the continuous stochastic process up to time t, { X (s), 0 ≤ s < t}, the conditional probability density Pn (t) ≡ E [P[ N (t) = n|{ X (s), 0 ≤ s < t}]] satisfies the master equation dPn = λ( X (t))[ Pn−1 (t) − Pn (t)], dt (4.36) which has the solution Pn (t) = with Λ(t) = (Λ(t))n −Λ(t) e , n! Z t 0 λ( X (t0 ))dt0 . (4.37) (4.38) We now observe that the rate function Λ(t) is itself stochastic with respect to different realizations { X (s), 0 ≤ s < t}. Therefore, in order to determine the probability Pn (t) that 80 the number of events in [0, t) satisfies N (t) = n, it is necessary to average with respect to these different realizations. That is, Pn (t) = EX [P[ N (t) = n|{ X (s), 0 ≤ s < t}] Z t n Z t 1 λ( X (s))ds , λ( X (s))ds exp − = EX n! 0 0 (4.39) where EX denotes expectations with respect to the stochastic process X. Introducing the characteristic function h GΛ(t) (z) = EX e izΛ(t) i , (4.40) it immediately follows that Pn (t) is related to the nth derivative of GΛ(t) (z): Pn (t) = (−i )n (n) GΛ (i ). n! (4.41) Furthermore, we can express the characteristic function for N (t) in terms of GΛ(t) : n h i 1 iz izN (t) izn −Λ(t) e Λ(t) e GN (t) ( z ) ≡ E e = ∑ e Pn (t) = ∑ EX n! n ≥0 n ≥0 " # n h iz i 1 = EX ∑ eiz Λ(t) e−Λ(t) = EX ee Λ(t) e−Λ(t) n≥0 n! = GΛ(t) (i − ieiz ). (4.42) It immediately follows that dGN (t) (z) E[ N (t)] ≡ −i dz z =0 dGΛ(t) (i − ieiz ) = −i dz = EX [Λ(t)]. (4.43) z =0 In order to determine more general statistics of the DSPP such as the covariance, we need to determine the joint characteristic function of a finite set of variables { N (t1 ), . . . , N (tm )}. This can be achieved using the notion of a characteristic functional [10]. The latter is defined according to Z Φ N [v] ≡ E exp i T 0 v(σ)dN (σ) , (4.44) for fixed T, where v is a real-valued function and the integral is a counting integral, that is, Z T 0 N (T ) v(σ)dN (σ) = ∑ v ( ωi ), (4.45) i =1 with ωi denoting the occurrence times of the DSPP. Expectation is taken with respect to both stochastic processes N (t), X (t). In order to evaluate the characteristic functional, 81 we first condition on a particular realization { x (t), 0 ≤ t ≤ T } of the stochastic process X (t) over the time interval [0, T ]. We write the corresponding conditioned characteristic functional as ( Φ N [v| x ] ∼ E exp i M ∑ v(σk )∆N (σk ) k =1 !) . (4.46) Following along analogous lines to the analysis of path-integrals [12], we have discretized time into M intervals of size ∆σ, and expectation is taken with respect to the inhomogeneous Poisson process with intensity λ(t) = λ( X (t)). (We are assuming that the limit M → ∞, ∆σ → 0 with M∆σ = T is well-defined.) It follows that Φ N [v| x ] ∼ M ∏ [(1 − λ(σk )∆σ) + λ(σk )∆σ exp (v(σk ))] k =1 M ∼ ∏ exp k =1 h e iv(σk ) i − 1 λ(σk )∆σ ∼ exp M ∑ k =1 h e iv(σk ) i ! − 1 λ(σk )∆σ . (4.47) If we now retake the continuum limit , we see that Φ N [v| x ] = exp Z T 0 h e iv(σ ) i − 1 λ(σ)dσ . (4.48) Finally, taking expectation with respect to the stochastic process X (t) yields the result Φ N [v] = EX exp Z T 0 h e iv(σ) i − 1 λ(σ)dσ . (4.49) Now take v(σ) to be the following piecewise function [10]:  m 0 ≤ σ < t1 , ∑ i =1 α i ,    m  α , t  1 ≤ σ < t2 ,  ∑ i =2 i . .. .. v(σ) = .    t m −1 ≤ σ < t m ,  αm ,   0, tm ≤ σ < T, (4.50) where 0 < t1 < t2 < · · · < tm < T. From equation (4.44), the corresponding characteristic functional is Φ N [v] = E {exp[i (α1 + · · · + αm ) N (t1 ) + i (α2 + · · · + αm )( N (t2 ) − N (t1 )) + · · · + iαm ( N (tm ) − N (tm−1 ))]} = E {exp[i (α1 N (t1 ) + · · · + αm N (tm ))]} = GN (t1 ),...,N (tm ) (α1 , . . . , αm ), (4.51) 82 where GN (t1 ),...,N (tm ) is the joint characteristic function of ( N (t1 ), . . . , N (tm )). On the other hand, from equation (4.49), we have Φ N [v] = EX i (α1 +···+αm ) t1 iαm Z tm λ(σ)dσ + · · · + e − 1 λ(σ)dσ t m −1 n h io = EX exp ei(α1 +···+αm ) − ei(α2 +···+αm ) Λ(t1 ) + · · · + eiαm − 1 Λ(tm ) = GΛ(t1 ),...,Λ(tm ) −iei(α1 +···+αm ) + iei(α2 +···+αm ) , . . . , i − ieiαm , (4.52) exp e −1 Z 0 where GΛ(t1 ),...,Λ(tm ) is the joint characteristic function of (Λ(t1 ), . . . , Λ(tm )) For the sake of illustration, consider the case m = 2 and the covariance function R N (t1 , t2 ) = E[ N (t1 ) N (t2 )] − E[ N (t1 )]E[ N (t2 )] ∂2 GN (t1 ),N (t2 ) (α1 , α2 ) − EX [Λ(t1 )]EX [Λ(t2 )] =− ∂α1 ∂α2 α1 = α2 =0 ∂2 GΛ(t1 ),Λ(t2 ) (−iei(α1 +α2 ) + ieiα2 , i − ieiα2 ) =− − EX [Λ(t1 )]EX [Λ(t2 )] ∂α1 ∂α2 α1 = α2 =0 = EX [Λ(t1 )Λ(t2 )] + EX [Λ(t1 )] − EX [Λ(t1 )]EX [Λ(t2 )] = RΛ (t1 , t2 ) + EX [Λ(t1 )]. (4.53) Expressing Λ(t) in terms of the intensity finally gives R N ( t1 , t2 ) = Z t1 Z t2 0 0 0 0 Rλ (τ, τ )dτdτ + Z t1 0 EX [λ(τ )]dτ, t1 < t2 , (4.54) where Rλ (τ, τ 0 ) = EX [λ(τ )λ(τ 0 )] − EX [λ(τ )]EX [λ(τ 0 )]. 4.2.2 (4.55) Calculation of Mean and Variance of IFT Numbers within Flagellum The above analysis shows that determining the first-order and second-order statistics of the number N (t) of injected IFT particles requires calculating the corresponding statistics of the stochastic intensity λ( X (t)) = η MX (t), where X (t) is the fraction of bound binding sites in the basal body. Thus, calculating the mean and covariance of the intensity reduces to determining these quantities for the solution of the SDE (4.26). Performing the Ito integral shows that X (t) = X ∗ h 1−e −Γt i + X0 e −Γt + r Γ M Z t 0 0 e−Γ(t−t ) q X ∗ + ΘX (t0 )dW (t0 ), (4.56) 83 for a fixed initial condition X (0) = X0 with X ∗ given by equation (4.3) and k− − k+ B . k− + k+ B (4.57) h i h X (t)i = X ∗ 1 − e−Γt + X0 e−Γt , (4.58) Γ = k+ B + k− , It follows that and so E[ N (t)] = ηM For t1 < t2 , we also have Z t 0 Θ= i X0 − X ∗ h −Γt , h X (τ )idτ = ηM X t + 1−e Γ ∗ (4.59) R X (t1 , t2 ) ≡ h[ X (t1 ) − h X (t1 )i][ X (t2 ) − h X (t2 )i]i q Z Z Γ t1 − Γ ( t1 − s ) t2 − Γ ( t2 − s 0 ) 0 ∗ ∗ 0 = [ X + ΘX (s)][ X + ΘX (s )]dW (s)dW (s ) e e M 0 0 Z Γ t1 −Γ(t1 +t2 −2s) ∗ e [ X + Θh X (s)i]ds = M 0 Z t1 Z t1 Γ − Γ ( t1 + t2 ) ∗ 2Γs ∗ Γs = e X (1 + Θ ) e ds + Θ[ X0 − X ] e ds M 0 0 Θ[ X − X ∗ ] X ∗ (1 + Θ ) − Γ [ t2 − t1 ] 0 e − e − Γ [ t1 + t2 ] + e−Γt2 − e−Γ[t1 +t2 ] . (4.60) = 2M M The case t1 > t2 is taken into account by exchanging t1 , t2 . Hence, setting Rλ = (η M)2 R X in equation (4.54) gives R N ( t1 , t2 ) = Mη 2 X ∗ (1 + Θ) 2 Z t1 Z t2 0 0 0 0 (e−Γ|τ −τ | − e−Γ[τ +τ ] )dτdτ 0 Z t1 Z t2 0 0 (e−Γ max{τ,τ } − e−Γ[τ +τ ] )dτdτ 0 + E[ N (t1 )] + Mη 2 Θ[ X0 − X ∗ ] 0 0 ∗ 2 X (1 + Θ ) ∗ = Mη A(t1 , t2 ) + Θ[ X0 − X ]B(t1 , t2 ) + E[ N (t1 )], (4.61) 2 where A(t1 , t2 ) = i i 2t1 2 h 1 h − 2 1 − e−Γt1 + 2 1 − e−Γ(t2 −t1 ) − e−Γt1 + e−Γt2 Γ Γ Γ ih i 1 h − 2 1 − e−Γt1 1 − e−Γt2 , Γ (4.62) and B(t1 , t2 ) = i t h i ih i 2 h 1 h 1 −Γt1 −Γt1 −Γt2 −Γt1 −Γt2 1 − e − e + e − 1 − e 1 − e . Γ2 Γ Γ2 In particular, setting t1 = t2 = t yields the variance ∗ X∗ 2 X A(t) + Θ[A(t) − 2B(t)] + ΘX0 B(t) + E[ N (t)], Var[ N (t)] = Mη 2 2 (4.63) (4.64) 84 with A(t) = i i2 2 h 1 h 2t − 2 1 − e−Γt − 2 1 − e−Γt Γ Γ Γ (4.65) and B(t) = i 2t i2 2 h 1 h −Γt −Γt −Γt . 1 − e − e − 1 − e Γ2 Γ Γ2 (4.66) Note that A(t) > 2B(t) > 0 for all t > 0 and Var[ N (t)] > E[ N (t)]. (4.67) The latter is a basic property of DSPPs, namely, that the variance is greater than a Poisson process with intensity given by the mean of the stochastic intensity - a feature known as over-dispersion. Ignoring transient statistics, we conclude that for large t, E[ N (t)] ∼ λ0 t, λ0 = η MX ∗ , (4.68) and Var[ N (t)] ∼ λ1 t, λ1 = λ0 + η 2 MX ∗ 1+Θ . Γ (4.69) Using a similar analysis to section 4.1.1, it follows that E[ F (t)] ∼ λ0 T, 4.2.3 Var[ F (t)] ∼ λ1 T. (4.70) Numerical Results We simulate the DSPP using a thinning algorithm [75] as follows. Consider a nonhomogeneous Poisson process on the time interval [0, T ] with rate function λ(t), and assume there exists a constant λ∗ such that λ∗ ≥ λ(t) on [0, T ]. To simulate the nonhomogeneous Poisson process, first consider the homogeneous Poisson process with rate λ∗ . We now generate a sequence of times T1 , T2 , ..., Tm , for m ∈ N with 0 < T1 < T2 < ... < Tm ≤ T, with Ti , i = 1, ..., m corresponding to the time of the ith injection of IFTs docked at the basal body into the flagellum. To obtain the sequence of injection times for the nonhomogeneous Poisson process with rate λ(t), we accept each Ti generated from the homogeneous Poisson process with probability λ( Ti )/λ∗ . The resulting sequence of injection times corresponds to the nonhomogeneous Poisson process with rate function λ(t). For a rigorous proof of this, see [75]. For our particular model, we employ the thinning algorithm by utilizing the following procedure: 85 1. Generate a stochastic trajectory X (t) according to equation (4.26) on the interval [0, T ]. 2. Compute λ( X (t)) = ηMX (t) and let λ∗ = max(λ( X (t))). 3. Generate a sequence of times T1 , T2 , ..., Tm with 0 < T1 < T2 < ... < Tm ≤ T from an exponential distribution with parameter λ∗ . 4. For each Ti , i = 1, .., m, generate a random number Ui distributed uniformly on the interval [0, 1]. 5. If λ( X ( Ti )) λ∗ ≥ Ui accept Ti as a firing time generated by the nonhomogeneous Poisson process with rate λ( X (t)). Otherwise, do not include Ti as a firing time generated by the nonhomogeneous Poisson process. In Figure 4.5 we show histograms for the number of IFTs injected into the flagellum on the time interval [0, T ] generated from the numerical procedure elucidated above for cases where M = 200 and M = 20. Both of these stochastic processes depict over-dispersion. We can see that the histogram in the small M case is considerably noisier than the histogram in the large M case, which is relatively smooth. This dichotomy in the histograms depicts the beginning of the breakdown of the system-size expansion in the small M regime. In p Figure 4.6, we show plots of h N ( T )i and the coefficient of variation, Var[ N ( T )]/h N ( T )i, 0.06 0.2 a.) b.) 0.18 0.05 0.16 0.14 0.04 0.12 0.03 0.1 0.08 0.02 0.06 0.04 0.01 0.02 0 20 30 40 50 N(t) 60 70 0 350 400 450 500 N(t) Figure 4.5. Histogram depicting N ( T ) for L = 10µm over 1000 trials. (a) M = 20. Here, h N ( T )i = 43.477 and Var[ N (5)] = 50.599. (b) M = 200. Here, h N ( T )i = 433.284 and Var[ N ( T )] = 435.989. Both histograms depict over-dispersion. Other parameter values are B = 100, σ = κ = 5, η = 1, γ = 20, k − = 1, τ = 1, v = 2, k + = 1. 500 0.07 450 0.065 400 0.06 CV Total IFT 86 350 0.055 300 0.05 250 0.045 200 a.) 0 5 10 Ciliary Length [µm] 15 0.04 0 5 10 Ciliary Length [µm] b.) 15 Figure 4.6. Simulation of the DSPP with using equation (4.26) averaged over 100 trials. (a) Plot of h N ( T )i. Simulations with error bars shown in blue. Analytical curve shown in red. (b) Plot of coefficient of variation versus ciliary length. Blue points are simulation results. Analytical curve shown in red. Other parameter values are η = 1, σ = κ = 5, γ = 20, B = 100, M = 200, k − = 10,k + = 1, τ = 1, v = 2. versus ciliary length, with stochastic trajectories generated by the aforementioned numerical procedure, and compare them to the analytical results given in equations (4.59) and (4.64) (In particular, to compute the coefficient of variation, we take the square root of equation (4.64) and divide by equation (4.59).) The numerical results are in very good agreement with the analytical curves. 4.3 Modified Stochastic Model for Small M Our formulation of flagellar length control as a DSPP relied on the assumption that the number M of binding sites within the basal body is sufficiently large so that the binding process is independent of the Poisson process (but not vice versa). On the other hand, when M = O(1), we have to consider the joint stochastic process that simultaneously keeps track of the number of bound IFTs N (t) and the number M(t) particles injected into the flagellum. The master equation for Pm,n (t) = P( M(t) = m, N (t) = n| M (0) = m0 , N (t) = 0) is dPm,n = k + B( M − m + 1) Pm−1,n (t) + k − (m + 1) Pm+1,n (t) dt − [k + B( M − m) + k − m] Pm,n (t) + η (m + 1) Pm+1,n−1 (t) − ηmPm,n (t). The mean number of particles injected in the interval [0, t] is (4.71) 87 M h N (t)i = ∑ ∑ nPm,n (t), (4.72) m =0 n ≥0 and the mean number of particles within the flagellum at time t is h F (t)i = h N (t)i − h N (t − T )i. (4.73) Unfortunately, is not possible to obtain exact analytical solutions to the full master equation, so we will investigate the effects of small M using computer simulations. However, before presenting our numerical results, it is important to understand how the master equation (4.71) is related to the DSPP master equation (4.36) for large M. The first step is to carry out a partial system-size expansion of equation (4.71) with respect to the number of bound IFTs m along similar lines to section 4.2.2. Setting Pn ( x, t) = Pm,n (t) with m = Mx, we obtain the differential Chapman-Kolmogorov (CK) equation 1 ∂2 ∂ ∂Pn ( x, t) = − [ A( x ) Pn ( x, t)] + [ B( x ) Pn ( x, t)] ∂t ∂x 2M ∂x2 + λ( x ) Pn−1 ( x, t) − λ( x ) Pn ( x ), (4.74) with λ( x ) = η Mx. The next step is to show that Pn ( x, t) = E[ Pn (t)1X (t)= x ], (4.75) where Pn (t) is the solution to the stochastic master equation (4.36) and the expected value is taken with respect to the Brownian process X (s), 0 ≤ s < t. We will establish this using the theory of Brownian functionals. 4.3.1 Feynman-Kac Formula and Reduction to a DSPP Let X (t) ∈ R represent Brownian motion, such as the solution to the SDE (4.26). A Brownian functional over a fixed time interval [0, t] is formally defined as a random variable Λ given by Λ= Z t 0 U ( X (τ ))dτ, (4.76) where U ( x ) is some prescribed function or distribution such that Λ has positive support. In our case, we take U ( X (t)) = λ( X (t)) = η MX (t) and identify Λ(t) with the rate function of the DSPP. Since X (t), t ≥ 0 is a Wiener process, it follows that each realization of a Brownian path will typically yield a different value of Λ, which means that Λ will be 88 distributed according to some probability density P(Λ, x, t| x0 , 0) for X (0) = x0 and X (t) = x. The statistical properties of a Brownian functional can be analyzed using path integrals, and lead to the well-known Feynman-Kac formula [54]. For a general review of Brownian functionals and their applications, see reference [79]. Here we briefly indicate the steps in the derivation of the Feynman-Kac formula. Motivated by the analysis of DSPPs in section 4.2.1, introduce the characteristic functional G (z, x, t| x0 , 0) = Z ∞ 0 e(iz−1)Λ P(Λ, x, t| x0 , 0)dΛ. (4.77) Using the classical path-integral representation of Brownian motion, we have Z x(t)= x Z ∞ Z t G (z, x, t| x0 , 0) = e(iz−1)Λ δ Λ− λ( x (τ ))dτ P[ x ]D[ x ] dΛ 0 x (0)= x0 0 Z x(t)= x Z t = exp (iz − 1) λ( x (τ ))dτ P[ x ]D[ x ] x (0)= x0 0 = E exp (iz − 1) Z t 0 λ( x (τ ))dτ x(t)= x . (4.78) x(t+∆t)= x (4.79) x (0)= x0 Now note that G (z, x, t + ∆t| x0 , 0) = E exp (iz − 1) ≈ Z Z t+∆t 0 λ( x (τ ))dτ p(∆W )E exp (iz − 1) x (0)= x0 Z t 0 λ( x (τ ))dτ × exp ((iz − 1)λ( x )∆t) d∆W. x(t)= x−∆x (4.80) x (0)= x0 (4.81) We have split the time interval [0, t + ∆t] into two parts [0, t] and [t, t + ∆t] and introduced the intermediate state x (t) = x − ∆x with; see equation (4.26), r B( x − ∆x ) ∆W (t). ∆x = A( x − ∆x )∆t + M (4.82) It follows that G (z, x, t + ∆t| x0 , 0) = e(iz−1)λ(x)∆t Z p(∆W ) G (s, x − ∆x, t| x0 , 0)d∆W ∂ h∆x i G (z, x, t| x0 , 0) ∂x ∂2 2 + 2 h∆x i G (z, x, t| x0 , 0) + . . . . ∂x = e(iz−1)λ(x)∆t G (z, x, t| x0 , 0) + (4.83) (4.84) (4.85) Using the fact that for the SDE (4.26) h∆x i = A ( x ), ∆t→0 ∆t lim h(∆x )2 i = B ( x ), ∆t ∆t→0 lim (4.86) 89 we obtain, in the limit ∆t → 0, the desired Feynman-Kac formula in the form of a modified Fokker-Planck equation: ∂G ∂A( x ) G ∂2 B( x ) G =− + + (iz − 1)λ( x ) G. ∂t ∂x ∂x2 (4.87) Note that G satisfies the initial condition G ( x, t0 | x0 , t0 ) = δ( x − x0 ). (4.88) The final step in connecting the master equation (4.36) with the CK equation (4.74) is to note that 1 E[ Pn (t)1X (t)= x ] = n! ∂ −i ∂z n G (z, x, t| x0 , 0) z =0 ≡ Gn ( x, t). (4.89) Therefore, differentiating equation (4.87) with respect to z yields ∂Gn ∂A( x ) Gn ∂2 B( x ) Gn − λ( x ) Gn + λ( x ) Hn , =− + ∂t ∂x ∂x2 where 1 Hn ( x, t) = n! ∂ −i ∂z n {izG (z, x, t| x0 , 0)} . (4.90) (4.91) z =0 It is straightforward to establish that Hn ( x ) = Gn−1 ( x ) and thus, equation (4.90) is equivalent to equation (4.74) on identifying Pn ( x, t) with Gn ( x, t). 4.3.2 Numerical Results For numerical simulation of the system in the small M regime, we use a continuoustime Monte Carlo algorithm based on the Gillespie algorithm [35]. Plots for h N ( T )i and the coefficient of variation versus ciliary length are shown in Figure 4.7. In the large M regime (Figure 4.7a and Figure 4.7b), we find that the stochastic curve generated from the chemical master equation (4.74) is in good agreement with the analytical results from equations (4.59) and (4.64). This is numerical evidence for the equivalency of the chemical master equation and the DSPP in the large M limit. Figure 4.7c and Figure 4.7d depict the mean and coefficient of variation of injected IFTs versus ciliary length in the small M regime. As noted previously, the system-size expansion breaks down in this regime, whereas the chemical master equation (4.74) holds for all M values. This is indicated by the difference between the results of the stochastic simulations and the analytical curves in Figure 4.7c and Figure 4.7d. 500 0.07 400 0.06 CV Total IFT Injected 90 300 200 0.05 a.) 0 b.) 5 10 0.04 15 5 10 15 5 10 15 0.4 12 10 CV Total IFT Injected 14 0 0.35 8 6 4 0.3 c.) 0 d.) 5 10 15 Ciliary Length [µm] 0 Ciliary Length [µm] Figure 4.7. Simulation results of equation (4.74) averaged over 200 trials. Figures (a) and (b) have M = 200. Figures (c) and (d) have M = 5. Other parameter values are the same as Figure 4.6. 4.4 Discussion In this chapter, we presented a description of flagellar length control based on a Doubly Stochastic Poisson Process (DSPP) wherein we assumed the injection times of IFT into a flagellum from the basal body are given by a nonhomogeneous Poisson process whose rate is based on the number of IFTs bound at the basal body. The number of IFTs bound to the basal body in turn evolves according to a stochastic birth death process, hence rendering the Poisson rate of injection times of IFT into the flagellum stochastic. In the case where the number of binding sites on the basal body is sufficiently large, we invoked the system-size expansion to derive a Langevin equation describing the time evolution of the density of occupied binding sites. We showed how the statistics of the number of IFTs injected into the flagellum are related to the density of occupied sites on the basal body in the large M limit using a characteristic functional argument. In particular, we showed how the solution of the DSPP exhibits over-dispersion. We then relaxed the condition that the number of 91 binding sites is large and described the full stochastic process using a chemical master equation. We then invoked the Feynman-Kac formula to describe how the chemical master equation reduces to the DSPP. In particular, we have shown that the underlying dynamics of our DSPP are completely described by a chemical master equation. We also point out that our doubly stochastic description of IFT injection into the flagellum supports the existence of a unique equilibrium flagellar length, in agreement with the model in [82]. We show this in Figure 4.8, where we plot h N ( T )i/L versus ciliary length. We find the curve monotonically decreases, indicating that for any constant ζ ≡ 2V/av (with the parameters defined in the introduction to this chapter), there is a unique length L∗ such that h N ( T )i/L∗ = ζ. There are a number of possible extensions of our work. In this model, we have assumed the motion of the IFTs along the flagellum is ballistic, whereas in reality, the motion is more random. We can model this by describing the motion of the IFTs by an advection-diffusion equation, for example, along the lines of Chapter 2. It would also be interesting to see the impact of allowing for the τ parameter, which reflects the amount of time spent at the tip of the flagellum by the IFT, to be a random variable. Generalizing more, another aspect to study would be to consider the cases where the number of particles available to bind to a <N(T)>/L 300 200 ζ 100 0 0 5 10 15 Ciliary Length [µm] Figure 4.8. Plot showing relationship between average number of IFTs in flagellum and ciliary length, and the existence of a unique stationary flagellum length for some constant ζ. Parameter values are the same as Figure 4.6. 92 basal body (in our model this is represented by B) is finite and is divided between two or more competing flagella [36, 94] . CHAPTER 5 AXONAL LENGTH SENSING The problem of length control highlighted early in Chapter 4 is particularly acute for the axons of neurons, which exhibit the most significant size differences of any cell type, ranging from several microns to a meter in humans. It is likely that different growth mechanisms operate at different stages of development [37]. For example, the initial growth of an axon is determined by preprogrammed transcription factor levels (quantal synthesis) [70], whereas the interstitial growth rates of axons that have connected to their targets is driven by stretching of the organism [112]. A major question is whether or not there is an intrinsic length sensor that can coordinate between transcriptional and metabolic process controlled by the nucleus and the differential growth and maintenance of axonal length. In vitro experimental studies of axonal growth in a variety of neuronal types support the existence of intrinsic length sensors [38, 55, 90, 106], but the underlying mechanisms are still largely unknown. Given the lengths involved, it is unlikely that a diffusion-based mechanism or a molecular ruler such as a microtubule can be involved. Recently, Rishal et al. [3, 102] have proposed a bidirectional motor transport mechanism for cellular-length sensing in axons, which would form the front-end of a length control mechanism that is distinct from [I]-[III] listed above. A schematic illustration of the motor-based model is shown in Figure 5.1. An anterograde signal is transported by kinesin motors from the cell body to the tip of the growing axon, where it activates the dynein-mediated transport of a retrograde signal back to the cell body. The retrograde signal then represses the anterograde signal via negative feedback, resulting in an oscillating retrograde signal whose frequency decreases with axon length. If axonal growth rates are correlated with this frequency, then spatial information regarding the length of the axon can be communicated to the cell body, where the frequency-dependent activation of transcription factors [21] could regulate axonal growth. One major prediction of the model 94 anterograde signal soma growth cone retrograde signal Figure 5.1. Schematic diagram of bidirectional motor-transport mechanism for axonal length sensing hypothesized by Rishal et al. [102]. A kinesin-based anterograde signal activates a dynein-based retrograde signal that itself represses the anterograde signal via negative feedback. The frequency of the resulting oscillatory retrograde signal decreases with axonal growth. based upon computer simulations is that reducing either anterograde or retrograde signals (by the partial knockdown of kinesin or dynein motor activity) should increase axonal length. This prediction has been confirmed experimentally in peripheral sensory neurons [102]. Note that a previous model of Kam et al. [55] is inconsistent with the experimental data. The earlier model assumes that the unidirectional transport of a retrograde signal by dynein motors maintains axonal growth until the signal at the cell body becomes too weak due to a constant rate of signal loss en route. In this case, the partial knockdown of motor activity would lead to shorter axons. The experimental results therefore provide circumstantial evidence for frequency-encoded axonal length. Such cellular behavior has been shown to exist in the context of protein production in response to the gonadotropic releasing hormone (GnRH), which pulses at various frequencies over time [65]. Distinct frequencies have been observed to induce the production of disparate proteins. This phenomenon was mathematically analyzed in reference [65]. The results suggest that cellular decoding of frequency-encoded information is possible due to the difference in time scales for gene activity and protein lifetime. Even more interestingly, it has been shown that cells are able to keep protein levels with less variability in response to a pulsatile signal as 95 opposed to a constant signal [117]. In this chapter, we develop a mathematical model of the biophysical mechanism proposed by Rishal et al. [102], in order to carry out a more systematic investigation of the dynamical process generating oscillations, and how the oscillation frequency depends on various biophysical parameters. We first consider a simple delay-differential equation with negative feedback that models the chemical signals at the somatic and distal ends of the axon (see section 5.1.1). The molecular motors are not modeled explicitly; rather their active transport is assumed to introduce a discrete delay that varies linearly with axonal length. We show how oscillations arise at a critical axonal length via a Hopf bifurcation, and obtain a length-dependent frequency consistent with the previous computational model (see section 5.1.2). We then construct a system of advection-diffusion equations that couple the chemical signaling with the active transport of kinesin and dynein motors (see section 5.1.3). Each advection-diffusion equation is an effective mean field equation for the transport of a population of motors of a given type, which randomly switch between a motile state (bound to a microtubule) and a diffusive state (unbound to a microtubule) (see section 1.1). In section 5.2, we use Green's functions to show how the advection-diffusion model is structurally similar to the delayed-feedback model. Now, however, the discrete delay is replaced by a distribution of delays. In section 5.3, we use numerical simulations to confirm that the PDE model supports a similar length-dependent frequency to the delayed feedback model. We also show how knockdown of either motor type increases the frequency, thus leading to longer axons as found experimentally. One prediction of our model is that the critical axonal length at which oscillations first occur increases with the diffusivity D, but the frequency of oscillations beyond criticality is relatively insensitive to D. Although the diffusion term partially captures the stochastic nature of motor transport, there are additional levels of stochasticity that are not captured by the mean-field model. By carrying out numerical simulations of (i) the computational model of Rishal et al. [102] and (ii) a stochastic version of our advection-diffusion model, we show that there are fluctuations in the frequency of the oscillatory signal whose coefficient of variation (standard deviation/mean) increases monotonically with length. This suggests that the proposed mechanism for axonal length sensing could break down for long axons. In sections 5.4 and 5.5, we explore the decoding of the oscillatory chemical signal. First 96 we imagine feeding the oscillating retrograde signal from this model into a feed-forward serial gene network. We find that the mean protein output from such a gene network is a monotonically decreasing function of axonal length, thus providing a faithful representation of axonal length. If the protein output were thresholded, then this could provide a mechanism for axonal length control. We then investigate the impact of intrinsic noise due to finite copy numbers within the gene network on the relationship between mean protein output and axonal length, and obtain some error estimates for the critical axonal length. 5.1 Model Formulation We now elucidate our mathematical model of the biophysical mechanism proposed by Rishal et al. [102]. 5.1.1 Delayed Feedback Model Consider an axon of length L with x = 0 corresponding to the proximal end (adjacent to the cell body or soma) and x = L corresponding to the distal end (axonal tip). In this section, we will ignore the dynamics of L by exploiting the fact that axonal growth occurs much more slowly than the time-scales of motor transport. Under this adiabatic approximation , we can treat L as fixed and investigate the occurrence of oscillations in chemical signaling for fixed length. This will then be used to determine how the frequency of oscillations varies as a function of length. Let u E (t) denote the anterograde chemical signal at x = L and time t, which is transported by kinesin motors from the proximal end at x = 0. Similarly, let u I (t) denote the retrograde signal at x = 0 at time t, which is transported by dynein motors from the distal end x = L. For the moment, we will assume the simplest possible model of active transport, in which both types of motor travel at a constant speed v along the axon (via binding to polarized microtubules). This means that for given length L, there is delay τ = L/v between the production of a signal at one end and its arrival at the opposite end. This motivates the following delayed feedback model (see Figure 5.2): du E dt du I dt = I0 − γu E − WI f (u I (t − τ )), (5.1a) = −γu I + WE f (u E (t − τ )), (5.1b) 97 anterograde signal I0 WE axon UE WI UI x=0 x=L retrograde signal Figure 5.2. Schematic diagram of feedback model. See text for details. with γ a decay rate. For simplicity, we take the decay rate to be the same for both chemical signals. The weights WE and WI determine the strength of the positive and negative feedback terms based on some form of Michaelis-Menten kinetics so that f is given by the Hill function f (u) = un , K n + un (5.2) for dissociation constant K and Hill coefficient n. We take n = 4 and fix the scale of the weights WE , WI and the input I0 by setting K = 2. The constant input I0 determines the rate at which kinesin packets are released at x = 0 in the absence of negative feedback (WI = 0). In order to match up with the results of Rishal et al. [102], we take γ−1 = 100 sec. Since kinesin and dynein motor velocities v are on the order of 1 µm/sec, it follows that τγ = 1 corresponds to an axonal length of 100 µm. In the following, we fix the units of time by setting γ = 1. 5.1.2 Linear Stability Analysis In order to relate our model to the length-sensing mechanism hypothesized by Rishal et al. [102], we look for a periodic solution of the coupled system given by equation (5.1) and determine how the effective frequency ω of the solution (if it exists) depends on the delay τ and thus on the axonal length L. First, setting time derivatives to zero in equation (5.1) yields the steady-state solutions u∗E and u∗I : u∗E = I0 − WI f (u∗I ), u∗I = WE f (u∗E ). Linearizing equation (5.1) about the steady state yields the linear system (5.3) 98 y˙E = −y E − α I y I (t − τ ), (5.4a) = − y I + α E y E ( t − τ ), (5.4b) y˙I where, for P = E, I, y P (t) ≡ u P (t) − u∗P , y˙P ≡ dy P /dt, and α P ≡ WP f 0 (u∗P ). This has the solution y P (t) = eλt YP with λ determined from the eigenvalue equation (λ + 1)eλt YE = −α I YI eλ(t−τ ) , (5.5a) = α E YE eλ(t−τ ) . (5.5b) (λ + 1)eλt YI Following a standard analysis of delay differential equations, we determine necessary conditions for the emergence of a time-periodic solution via a Hopf bifurcation by setting λ = iω and YP = UP + iVP in equations (5.5a) and (5.5b). Equating real and imaginary parts in the resulting system yields a matrix equation A(UE , VE , U I , VI )> = 0 with   1 −ω α I cos (ωτ ) α I sin (ωτ )  ω 1 −α I sin (ωτ ) α I cos (ωτ )  . A=   −α E cos (ωτ ) −α E sin (ωτ ) 1 −ω α E sin (ωτ ) −α E cos (ωτ ) ω 1 (5.6) In order for (UE , VE , U I , VI )T to be nontrivial, we require the matrix A to have a zero √ determinant. It turns out that this holds if U I = VE = 0 and VI = ± α E /α I UE . We thus obtain the following conditions for a Hopf bifurcation: ω = cot (ωτ ), √ α E α I sin (ωτ ) = 1, (5.7) Clearly these conditions cannot be satisfied in the absence of a delay (τ = 0). Indeed, setting τ = 0 in equations (5.5a) and (5.5b) shows that there exists a pair of eigenvalues √ given by λ± = −1 ± −α I α E . Since the real part of λ± is always negative, it follows that the steady state (u∗E , u∗I ) is stable, and periodic solutions cannot exist in the absence of a √ delay. On the other hand, if α E α I > 1, then a pair of complex conjugate eigenvalues √ crosses the imaginary axis at a critical positive delay τc , which depends on α E α I ; see Figure 5.3a. Although this is not sufficient to guarantee the emergence of a stable periodic solution via a supercritical Hopf bifurcation for τ > τc , the existence of stable oscillations beyond the Hopf bifurcation point can be verified numerically as illustrated in Figure 5.4. Moreover, the frequency of the oscillation decreases monotonically with τ such that there is an approximately five-fold decrease in frequency when axonal length reaches ∼1000µm (see Figure 5.3b). This is in agreement with the computational model of Rishal et al. [102]. 99 10 (a) Frequency[10-4 Hz] Critical Delay (τc) 8 6 4 2 0 0 2 √αEαI 4 (b) 5 6 4 3 2 1 0 2 4 6 Length[100μm] 8 10 Figure 5.3. Plots verifying the existence of the relationship between axonal length and oscillation frequency. (a) Plot of critical delay τc as a function of the effective coupling √ parameter α E α I . Both are in units of 100 sec. (b) Frequency of periodic solutions plotted against axonal delay. Parameter values as in Figure 5.4. The Hopf Bifurcation condition given by equation (5.7) can only be satisfied if √ αE α I > 1. Since α P ≡ WP f 0 (u∗P ), it follows that the strengths WP of the chemical signals carried by kinesin and dynein, respectively, must be sufficiently strong and/or the Hill function must be sufficiently steep. The latter suggests that oscillations are facilitated if the interactions between the chemical signals and the opposing molecular motors are cooperative in nature, as determined by the value of the Hill coefficient n in equation (5.2). In conclusion, our simple mathematical model makes explicit the crucial role of negative feedback in the proposed frequency encoding mechanism for axonal length sensing, and provides an analytical framework for studying such a mechanism. However, as it stands, the model is too phenomenological. In particular, it does not explicitly take into account the motion of the molecular motors. In order to incorporate the latter, we now consider a spatially extended version of our model that takes the form of an advection-diffusion equation. 5.1.3 Advection-Diffusion Model Let c± ( x, t) denote the density of kinesin (+) and dynein (−) motors at position x along the track. A simple model of active motor transport is to assume the motor densities evolve according to an advection-diffusion equation [12, 98, 111]: ∂c2 ∂c+ ∂c+ = −v + D +2 , ∂t ∂x ∂x ∂c2− ∂c− ∂c− =v +D 2. ∂t ∂x ∂x (5.8a) (5.8b) 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 (a) E Signal I Signal 0 20 Chemical Signals Chemical Signals 100 40 60 80 Time[100 sec] 5 4 3 2 1 0 (b) 0 40 60 80 Time[100 sec] 20 6 (c) Chemical Signals Chemical Signals 6 100 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 100 (d) 5 4 3 2 1 0 20 40 60 Time[100 sec] 80 100 0 0 20 40 60 Time[100 sec] 80 100 Figure 5.4. Chemical signal oscillations in the delayed feedback model given by equation (5.1) for various values of the delay (in units of 100 sec): (a) τ = 0.2; (b) τ = 0.29; (c) τ = 0.75; (d) τ = 1.5. Other parameters values are n = 4, I0 = 6, WE = WI = 5.5 such that τc ≈ 0.25. These are supplemented by the following boundary conditions at the ends x = 0, L: D∂ x c+ ( L, t) = 0, D∂ x c− (0, t) = 0, (5.9) J+ (0, t) = JE (t), J− ( L, t) = J I (t), (5.10) and where we have introduced the fluxes J± ( x, t) = ±vc± ( x, t) − D ∂c± . ∂x (5.11) The boundary conditions given by equation (5.9) impose the condition that the Fickian contribution to the flux of motors exiting the axon is zero. We are also assuming that at the end x = 0, kinesin motors are injected in the anterograde direction at a rate JE (t), whereas at the end x = L, dynein motors are injected at a rate J I (t) in the retrograde direction. For 101 simplicity, we take the mean speed and diffusivity of both motor species to be the same. Although dynein motors tend to move more slowly than kinesin motors, they are of the same order of magnitude. Moreover, we obtain similar results if the differences between the motors are taken into account. Suppose that each kinesin motor complex carries an amount κ E of excitatory chemical signaling molecules XE and each dynein motor carries an amount κ I of inhibitory chemical signaling molecules X I . When kinesin motors reach x = L, they release the molecules XE , which then enhance the injection rate of mobile dynein motors, whereas when dynein motors reach x = 0, they release the molecules X I , which then reduce the injection rate of mobile kinesin motors. We thus take JE (t) = v( IE − w I f [u I (t)]), (5.12a) J I (t) = −vwE f [u E (t)], (5.12b) where u E (t) is the concentration of XE at x = L, u I (t) is the concentration of X I at x = 0, and IE is the flux of kinesin motors injected at the proximal end. The latter evolve according to the pair of equations du E = −γu E + κ E c+ ( L, t), dt du I = −γu I + κ I c− (0, t). dt (5.13a) (5.13b) In the case of pure ballistic motor transport at a fixed speed v and D = 0, the above model reduces to our delayed feedback model. First note that the solution to equations (5.8) have the simple form c+ ( x, t) = F+ (t − x/v), c− ( x, t) = F− (t + x/v), (5.14) with the functions F± determined by the boundary conditions equation (5.10) - the boundary conditions (5.9) are not needed for these quasilinear differential equations. Thus F+ (t) = IE − w I f [u I (t)], F− (t + L/v) = wE f ([u E (t)]. (5.15) Substituting this into equations (5.13a) and (5.13b), we recover our previous model given by equation (5.1) with I0 = κ E IE , and w P = κ P WP for P = E, I. The spatially extended model can be analyzed along similar lines to the simpler model using Green's functions (see section 5.2). 102 5.2 Analysis of PDE Model Using Green's Functions In the case where the diffusion coefficient is nonzero, the solution to equations (5.8a,b) subject to the boundary conditions (5.9) and (5.10) is given by: c+ ( x, t) = Z L 0 G+ (ξ, 0; x, t)ψ+ (ξ )dξ + v Z t G+ (0, σ; x, t)[ IE − w I f (u I (σ))]dσ, (5.16) Z t (5.17) 0 and c− ( x, t) = Z L 0 G− (ξ, 0; x, t)ψ− (ξ )dξ + v 0 G− ( L, σ; x, t)wE f (u E (σ))dσ, where ψ± are the initial conditions for c± and G± are the Green functions for the respective advection-diffusion operators (see Appendix): ∞ G+ (ξ, σ, x, t) = ∑ e(−(Dλ + n v2 4D )( t − σ )) v e(− 2D (ξ − x)) (5.18a) n =1 √ p p 2D λn ×(sin ( λn x ) + cos ( λn x ))(sin ( λn ξ ) v √ p 2D λn + cos ( λn ξ )), v p ∞ G− (ξ, σ, x, t) = ∑ e−(Dη + n v2 4D )( t − σ ) v e 2D (ξ − x) (5.18b) n =0 √ p p 2D λn ×(sin ( λn x ) + cos ( λn x ))(sin ( λn ξ ) v √ p 2D λn cos ( λn ξ )), + v p √ √ where λn is an eigenvalue of the operator satisfying 4Dv λn cot ( λn L) = 4D2 λn − v2 . Substituting equations (5.16) and (5.17) into equations (5.13a) and (5.13b), respectively, yields: du E dt du I dt Z t = −u E + vκ E G+ (0, σ; L, t)( IE − w I f [u I (σ)])dσ , −∞ Z t = −u I + vκ I G− ( L, σ; 0, t)wE f [u E (σ)]dσ . −∞ (5.19a) (5.19b) We have taken the lower time limit to be t = −∞ in order to eliminate the transient terms. Equations (5.19a) and (5.19b) have the steady-state solution u∗E = vκ E G + ( L)( IE − w I f [u∗I ]), (5.20a) u∗I = vκ I G − ( L)wE f [u∗E ], (5.20b) 103 with G± ( L) = Z ∞ −∞ G± ( L, 0; 0, t)dt. (5.21) Linearizing (5.19a) and (5.19b) about the steady states gives Z t y˙E = −y E − vα I = −y I + vα E y˙I −∞ Z t −∞ G+ (0, σ; L, t)y I dσ, (5.22a) G− ( L, σ; 0, t)y E dσ. (5.22b) Introducing the causal Green's functions G± (t) = G± (t) H (t), where H (t) is the Heaviside function, we can take the upper time limit in the convolution integrals to be t = ∞. Fourier transforming the resulting linearized system using the convolution theorem then yields (iω + 1)YE = −vα I Gb+ (ω )YI , (5.23a) (iω + 1)YI = vα E Gb− (ω )YE , where Gb± (ω ) = Z ∞ −∞ G± (t)e−iωt dt, YP (ω ) = Z ∞ −∞ (5.23b) y P (t)e−iωt dt. (5.24) Equations (5.23a,b) are identical in form to equations (5.5a) and (5.5b) for λ = iω under the replacement e−iωτ → vGb+ (ω ). Thus, we can derive conditions for the occurrence of a Hopf bifurcation along similar lines to the delay differential equation model. That is, we take ω to be real and set YP = UP + iVP for P = E, I: (iω + 1)(UE + iVE ) = −vα I (U I + iVI )Gb+ (ω ), with (iω + 1)(U I + iVI ) = vα E (UE + iVE )Gb− (ω ), Gb+ (ω ) = ∞ ∑ 1 v2 v2 4D v e(( Dλn + 4D )σ) e(− 2D (ξ − x)) iω + Dλn + √ p p p 2D λn ×(sin ( λn x ) cos ( λn x ))(sin ( λn ξ ) v √ p 2D λn + cos ( λn ξ )), v ∞ v2 v 1 )σ) 2D (( Dλn + 4D Gb+ (ω ) = ∑ e e (ξ − x ) v2 n=0 iω + Dλn + 4D √ p p p 2D λn ×(sin ( λn x ) + cos ( λn x ))(sin ( λn ξ ) v √ p 2D λn + cos ( λn ξ )). v n =0 (5.25a) (5.25b) (5.26a) (5.26b) 104 We then equate real and imaginary parts and solve the resulting 4 × 4 matrix equation for the vector (UE , VE , U I , VI )T . 5.3 Results As shown in section 5.2, the advection-diffusion model is structurally similar to the simple delayed feedback model, except that the discrete delay τ = L/v is replaced by a distribution of delays given by a corresponding Green's function that depends on the axonal length L. This suggests that the advection-diffusion model will also exhibit oscillations beyond a critical length Lc , whose frequency decreases monotonically beyond Lc . This is indeed found to be the case, as illustrated in Figure 5.5. In contrast to the previous model, we can now also keep track of the density profile of the kinesin and dynein motors during a single cycle of the chemical oscillations. The variation in the density profiles at different points on the cycle is shown in Figure 5.6. Figure 5.6a shows a growing distribution of kinesin motors due to injection at the proximal end and negligible dynein. Subsequent excitation of the dynein motors by the chemical signal transported by the kinesin motors results in a growing dynein distribution (Figure 5.6b). This leads to inhibition of the kinesin motors (Figure 5.6c). The reduction in the kinesin motor density causes the density of dynein motors to diminish throughout the axon, which then allows the kinesin density to grow again (Figure 5.6d). One immediate issue that arises is how the emergence of oscillations and the lengthdependent frequency depend on the diffusivity D. We find that for sufficiently long axons, the frequency ω is approximately independent of D. On the other hand, the critical length Lc for the emergence of oscillations does depend on D. This is illustrated in Figure 5.7, which shows the variation in frequency as a function of length for different diffusivities. In each case, the frequency is a monotonically decreasing function of L, consistent with the delayed feedback model. Increasing the diffusion coefficient D also increases Lc so that it can result in the disappearance of the oscillations, as shown in Figure 5.8. 5.3.1 Knockdown of Molecular Motor Activity One of the key predictions of the mechanism proposed by Rishal et. al [102] is that when the number of motors in a given axon are inhibited or knocked down, the axon's 105 90 70 80 Excitatory Inhibitory 50 70 Chemical Signals Chemical Signals 60 60 40 50 40 30 30 20 20 10 0 10 (a) 0 400 600 800 Time[100sec] 0 1000 (c) 0 100 90 80 70 60 50 40 30 20 10 0 200 400 600 800 Time[100sec] (b) 1000 (d) Chemical Signals Chemical Signals 100 90 80 70 60 50 40 30 20 10 0 0 200 200 400 600 800 Time[100sec] 1000 0 200 400 600 800 Time[100sec] 1000 Figure 5.5. Chemical signal oscillations in the advection-diffusion model equation (5.8) for various axonal lengths: (a) L = 100 µm; (b) L = 500µm; (c) L = 1000µm; (d) L = 2000µm. Other parameters values are I0 = 10, wE = w I = 9, γ = 1, κ E = κ I = 1, v = 1µm s−1 , D = 0.1µm2 s−1 . length should grow due to a resulting increase in the frequency of the retrograde signal. We can model the partial knockdown of kinesin motors in an axon by decreasing the background rate of kinesin flux I0 . Similarly, we can model dynein knockdown by decreasing wE , the strength of the kinesin mediated excitation of dynein activity at the distal end. In both cases, we find that the required axon length for the normalized frequency to decay to a particular value is greater than when kinesin or dynein motors are knocked down to a lesser extent; see Figure 5.9. 5.3.2 Effects of Noise Modeling active motor transport in terms of an advection-diffusion equation is a meanfield treatment of the underlying stochastic transport mechanism, in which motors randomly switch between a motile state (bound to a microtubule) and a stationary or slowly diffusing state (unbound to a microtubule). Although the advection-diffusion equation 106 1 (b) 0 10 9 8 7 6 5 4 3 2 1 0 0.2 0.4 0.6 0.8 Position Along Axon [100μm] 1 (d) Motor Density (c) Motor Density 10 9 8 7 6 5 4 3 2 1 0 0 10 9 8 7 6 5 4 3 2 1 0 Motor Density Motor Density 10 9 (a) Dynein 8 Kinesin 7 6 5 4 3 2 1 0 0 0.2 0.4 0.6 0.8 Position Along Axon [100μm] 0.2 0.4 0.6 0.8 1 Position Along Axon [100μm] 0 0.2 0.4 0.6 0.8 Position Along Axon [100μm] 1 Figure 5.6. Spatial profiles for kinesin and dynein motors at different times during one cycle of period T = 2π/ω after transients have disappeared: (a) t = 0; (b) t = T/4; (c) t = T/2; (d) t = 3T/4. The initial condition for the kinesin motors is a hyperbolic secant function, whereas the initial condition for the dynein motors is zero. Here L = 100µm and other parameters are as in Figure 5.5. partially captures the stochastic nature of motor transport at the population level, there are additional sources of stochasticity not taken into account by the mean field model. First, the stochastic transport of an individual molecular motor is more accurately described by a differential Chapman-Kolmogorov (CK) equation, which determines the probability density that the motor is at a particular location and in a particular internal state (mobile or stationary) at a time t [12]. If the transition rates between the internal states are sufficiently fast compared to the hopping rate of motors along the filament, then an adiabatic reduction can be carried out to reduce the CK equation to a Fokker-Planck (FP) equation [89]. Furthermore, if the number of motors is sufficiently large and they move independently (i.e., no exclusion effects), then the concentration of motors can be represented by an advection-diffusion equation, which is obtained by multiplying the FP equation by the number of motors. It follows that additional stochastic effects arise in the case of a finite 6 D = 1μm2/s D = 0.1μm2/s 5 4 3 2 1 0 0 (a) 2 4 6 8 Axon Length [100μm] 10 Frequency[10-4 Hz] Frequency[10-4 Hz] 107 5 4 3 2 1 0 (b) 0 0.2 0.4 0.6 0.8 Axon Length [100μm] 1 Figure 5.7. Variation in frequency as a function of length for different diffusivities in the advection-diffusion model. Other parameters are the same as Figure 5.5. (a) Plot showing that the frequency is insensitive to D except at small axonal lengths. (b) Plot showing that the main affect of diffusivity is to modify the critical length at which a Hopf bifurcation occurs. (A zero frequency indicates that the system is operating below the Hopf bifurcation point so there are no oscillations.) number of motors and slow transition rates. Yet another possible source of noise is the random loss of chemical signals carried by the motors, or the failure of a motor to rebind to the track; these would also lead to additional decay terms in the mean-field model. The presence of various sources of noise suggests that the encoding of axonal length in terms of the frequency of a retrograde chemical signal could break down for long axons due to the accumulation of fluctuations. We will demonstrate this by presenting results of simulations of (i) a slightly modified version of the computational model of Rishal et al. [102] and (ii) a stochastic version of our advection-diffusion model. The former explicitly takes into account finite-size effects by tracking the motion of individual motors. However, rather than explicitly modeling the stochastic "stop-and-go" motion of molecular motors, Rishal et al. [102] assume that each motor undergoes ballistic motion with a constant velocity that is generated from an experimentally determined velocity distribution, one for kinesin and the other for dynein. We implemented their computational model using the same velocity distributions, but with a modified scheme for injecting motors at each end. That is, rather than injecting motors as packets of fixed size, we injected individual motors at an instantaneous rate given by JE (t) and J I (t), respectively; see equations (5.12a) and (5.12b). The background flux of kinesin motors at x = 0, JE , was taken to be 100 80 70 60 50 40 30 20 10 0 Chemical Signals Chemical Signals 108 (a) 0 100 200 300 400 Time[100sec] 500 80 70 60 50 40 30 20 10 0 Excitatory Inhibitory (b) 0 100 200 300 400 Time[100sec] 500 Figure 5.8. Chemical signals in the advection-diffusion model for (a) D = 0.1µm2 s−1 and (b) D = 100µm2 s−1 . Here, L = 500µm and other parameter values are as in Figure 5.5. motors per minute, each motor was assumed to carry one unit of chemical signal, and wE = w I = 9. The results of computer simulations are shown in Figure 5.10 averaged over 100 trials (see below). The mean oscillation frequency decreases monotonically with axonal length L, as in our advection-diffusion model, but the relative size of fluctuations increases with L. This is established by plotting the coefficient of variation (CV), which is the standard deviation over the mean, as a function of length. It can be seen that axons of length L = 1000µm have a CV≈ 0.12, indicating nontrivial noise levels. It is important to note that our advection-diffusion model is not a mean-field version of the above computational model, but is based on a more realistic model of the stop-and-go transport of molecular motors. (The random switching between motile and nonmotile states would generate the velocity distributions used in the Rishal et al. model.) Inclusion of finite-size effects is expected to generate some form of multiplicative noise terms in the underlying advection-diffusion model. We hope to explore this issue in more detail elsewhere. Here, we take a more brute-force approach by adding a constant white noise term to the right-hand side of equations (5.8): ∂c+ ∂t ∂c− ∂t ∂c2 ∂c+ + D +2 + µξ + ( x, t), ∂x ∂x ∂c2− ∂c− = v + D 2 + µξ − ( x, t), ∂x ∂x = −v with µ the noise strength, hξ ± ( x, t)i = 0 and (5.27a) (5.27b) 109 6 5 4 δ = 1.5 δ = 1.0 δ = 0.5 δ = 0.1 3 2 1 0 0 (a) 10 2 8 4 6 Length[100μm] 3.0 Frequency[10-4 Hz] Frequency[10-4 Hz] 7 wE = 16 wE = 8 wE = 4 2.5 2.0 1.5 1.0 0.5 0 0 (b) 10 2 8 4 6 Length[100μm] Figure 5.9. Variation in frequency as a function of length for (a) decreasing flux δ ≡ I0 − w I (representing knockdown of kinesin), and (b) decreasing excitatory coupling wE (representing knockdown of dynein). Other parameters are the same as Figure 5.5. hξ ± ( x, t)ξ ± ( x 0 , t0 )i = δ(t − t0 )δ( x − x 0 ), hξ ± ( x, t)ξ ∓ ( x 0 , t0 )i = 0. (5.28) The boundary conditions for this model are the same as in the deterministic advection -diffusion model. It is important to note that the above stochastic advection-diffusion model does not exactly conserve the number of molecular motors. However, conservation RL p is approximately satisfied, since 0 ξ ± ( x, t)dx ≈ 0. The units of µ are 1/ Length × Time. We find that the variance in the frequency of oscillations is approximately independent of length - the variance for different values of the noise strength µ is plotted in Figure 5.11a. Since the mean frequency is a monotonically decreasing function of length, it follows that the relative size of frequency fluctuations increases with length. This is illustrated in Figure 5.11b, where we plot the CV as a function of length for µ = 1. As in the computational model, the size of fluctuations increases monotonically with axonal length. Note that for µ = 1, the CV for L = 1000µm is about ten times larger than that obtained in the computational model; see Figure 5.10b. The two different models yield comparable-sized fluctuations for µ = 0.1. Irrespective of this, a major prediction of both stochastic models is that the variance in the frequency of oscillations grows with axonal length, indicating a significant degradation in the reliability of frequency-coding as a length-sensing mechanism for very long axons. Note that we use a heuristic method for measuring frequency fluctuations in Figure 5.10 110 -4 3.5 ×10 0.12 Coefficient of Variation (a) Frequency [Hz] 3 2.5 2 1.5 0.08 0.06 0.04 0.02 1 0.5 (b) 0.1 0 200 400 600 800 Axon Length [μm] 1000 0 200 400 600 800 Axon Length [μm] 1000 Figure 5.10. Effects of noise on frequency-encoding mechanism for axonal length sensing in the computational model of Rishal et al. [102] with modified end conditions (see text for details). (a) Average frequency versus axonal length. (b) Variance in frequency versus axonal length. The plots were obtained by running the simulations over 100 trials. and Figure 5.11. That is, we consider N trials for a fixed axonal length, and for each run, we look at the power spectrum of the retrograde signal. The frequency is assumed to be encoded by the peak of the power spectrum for ω > 0. The mean and variance of the set of peak frequencies across the N trials are then calculated for each axonal length. One limitation of this method is that identification of the peak of the spectrum is difficult at large axonal lengths. Nevertheless, a high CV is consistent with a broad power spectrum, as illustrated in Figure 5.12 in the case of the stochastic advection-diffusion equation. For short axons (L = 100µm) the spectrum is characterized by a sharp peak around the mean frequency of oscillations. On the other hand, the spectrum is much broader for long axons (L = 1000µm) so that it is difficult to extract the mean frequency. Of course, the length-scale at which noise becomes significant will depend on the value of the noise strength µ. However, the general trend is clear: any stochasticity in the motion of the molecular motors will lead to the accumulation of errors in the length-sensing mechanism as the length of the axon increases. 5.3.3 Numerical Methods and Parameter Values We simulated the advection diffusion model given by equation (5.8) using a Backward Euler time discretization. We used an upwind scheme for the spatial discretization associated with the advection term and a central difference scheme for the spatial discretization 111 Coefficient of Variation Var[Frequency] -11 5 ×10 4 3 2 1 -4 (a) -3 -2 log10(μ) -1 1 0.8 0.6 0.4 0.2 0 0 (b) 8 4 6 Length[100μm] 2 10 Figure 5.11. Effects of noise on frequency-encoding mechanism for axonal length sensing in the stochastic advection-diffusion model. (a) Variance in oscillation frequency (determined by the peak of the power spectrum) versus noise strength µ. (b) Variance in frequency versus axonal length. The plots were obtained by solving the stochastic partial differential equations over 1000 trials with µ = 1. of the diffusion term. Let Ujn be the numerical approximation to the true solution of c+ in equation (5.8) at the jth spatial lattice point and the nth time step, let Vjn be the numerical approximation to the true solution of c− in equation (5.8) at the jth spatial lattice point and the nth time step, let k be the time step, and let h be the spatial step. Then the finite difference scheme used was Ujn+1 − Ujn k Vjn+1 − Vjn k = −v =v Ujn+1 − Ujn−+11 h Vjn++11 − Vjn+1 h +D +D Ujn++11 − 2Ujn+1 + Ujn−+11 h2 Vjn++11 − 2Vjn+1 + Vjn−+11 h2 , , (5.29a) (5.29b) where v and D are the motor velocities and diffusion coefficients, respectively. For equations (5.19a) and (5.19b), we used an explicit scheme. To account for the Neumann and Robin boundary conditions, which are given by equations (5.9) and (5.10), respectively, we used ghost points. That is, if there are N spatial lattice points with Nh = L, we introduced points at the lattice point N + 1 and −1 to account for the boundary conditions: n n UN V1n − V−n 1 +1 − U N −1 = 0, = 0, 2h 2h U1 − U−1 vU0n − D = v( IE − w I f [u I (t)]), 2h V n − VNn −1 − vVNn − D N +1 = −vwE f [u E (t)]. 2h (5.30) (5.31) (5.32) 112 ×10 6 (a) 1.5 1 0.5 0 1 3 2 4 -4 Frequency [10 Hz] 5 10 Power/Frequency Power/Frequency 2 ×10 5 (b) 8 6 4 2 0 4 5 -5 Frequency [10 Hz] 2 1 3 6 Figure 5.12. Power spectra of the retrograde chemical signal generated by the stochastic advection-diffusion model of equation (5.27) for (a) L = 100µm and (b) L = 1000µm. Other parameters are as in Figure 5.5. We solve for the ghost lattice points in the schemes for the boundary conditions and substitute them into equation (5.29). The motor velocity v was chosen to be 1µm/sec based on generally known distributions of velocities for kinesin and dynein motors. Other parameters were chosen for suitable computation of equation (5.8). The one thing we made sure was to keep IE > w I to allow for the continued propagation of the solutions to equation (5.8). This corresponds to ensuring that the background flux of kinesin motors at the proximal end is sufficient to overcome the suppressive effects of the retrograde signal from the dynein motors. Choosing w I > IE causes an abrupt end to the solutions and is not realistic for our biological mechanism to function. Unless otherwise noted, parameter values were taken as follows: I0 = 10, wE = w I = 9, γE = γ I = 1, κ E = κ I = 1, v = 1µm s−1 , D = 0.1µm2 s−1 , L = 100µm, h = 0.1µm, k = 1sec. 5.4 Frequency Decoding by a Feedforward Gene Network Suppose that the oscillating retrograde signal from the delayed feedback model equations (5.1) reaches the nucleus of a given neuron and causes the rapid activation of some gene and subsequent production of some protein C fast fast u I (t) −→ Active Gene −→ C − → ∅, λ (5.33) 113 with λ decay rate. Note that this does not account for explicit kinesin and dynein dynamics. This motivates the following model for the dynamics of protein C, dc = h[u I (t)] − λc, dt (5.34) where c denotes the concentration of protein C and h[u] is a monotonically increasing function satisfying h → h∗ ∈ (0, ∞) as u → ∞. Define g(t) ≡ h[u(t)]. Then g is T-periodic, where T is the period of u I (t). Following [65], we obtain the time-dependent solution for c(t). Introduce the integrating factor eλt . Then, d λt (e c(t)) = g(t)eλt , dt Z t ⇒ c ( t ) = c ( t 0 ) e − λ ( t − t0 ) + t0 g(s)e−λ(t−s) ds. (5.35) We integrate over a period of u I (t) so that, for m ∈ N, c((m + 1) T ) = c(mT )e −λT = c(mT )e −λT + Z ( m +1) T +e mT −λT Z T g(s)e−λ((m+1)T −s) ds g(s)eλs ds. 0 (5.36) Equation (5.36) gives a recursive finite difference equation for c at integer multiples of the period of u I (t). For large m, we thus have c(mT ) = e−λT 1 − e−λT Z T g(s)eλs ds. 0 (5.37) Hence, c(t) converges to a T-periodic solution following any transient dynamics, as shown in Figure 5.13. In order to characterize the protein output in terms of the frequency ω of u I (t), we find the time average of c(t) post transience. This can be done by simply integrating equation (5.34) over a period of u I (t): c̄ = 1 λT Z T 0 g(s)ds ≡ ḡ . λ (5.38) Equation (5.38) is an intuitive result. It says that the average protein output from the feedforward serial network is equal to the ratio of the average protein activation rate to the protein decay rate. To make the relationship between c̄ and T more explicit, we perform the following. Assume that in the posttransient time regime, the maximum value of u I (t) is given by UM and that the minimum value is given by Um , and that the u I transitions from UM to P 114 Figure 5.13. Simulation of the feed forward serial network equation (5.34) in response to a retrograde signal from equation (5.1). (a) Retrograde signal being fed into gene network, τ = 5. (b) Convergence of the solutions of equation (5.34) to a T-periodic solution post transience. h[u] is taken to be the same function as f [u] defined in equation (5.2) multiplied by a factor of 1000, and we set λ = 0.01. Other parameter values are n = 4, I0 = 10, WE = WI = 9.5 such that τc ≈ 1.5. Um occur very quickly compared to other temporal dynamics. Further assume that h[u] is a Hill function with a large Hill coefficient, so that h[UM ] = A and that h[Um ] ≈ 0. Let η < T denote the amount of time for which u I (t) is at its maximum value in a given period, η = κT for 0 < κ < 1. Then, ḡ ≈ Aη/T and c̄ ≈ Aη . λT (5.39) Note that the assumptions made regarding u I (t) are consistent with the behavior of the retrograde signal for sufficiently long delays (see Figure 5.5). Equation (5.39) suggests that if η is constant, then the mean protein output c̄ is dependent on the frequency T −1 of the pulsatile retrograde signal u I (t), the protein decay rate λ, and the rate of protein activation A. In particular, c̄ is inversely related to T, and more specifically it is a monotonically decreasing function of T. In the context of the delayed feedback model, this means that c̄ is a monotonically decreasing function of axonal length L, as illustrated in Figure 5.14. Though the analytical representation of c̄ was obtained by making assumptions that simplified the analysis of equation (5.34), it is nevertheless a reasonable reflection of the true relationship between c̄ and L, which can be obtained numerically; see Figure 5.14. Note that if UM is sufficiently large, then A ≈ h∗ due to the saturating nature of h. What is more, changing the value of UM will not alter c̄ significantly 115 Mean Protein Output 9500 9000 8500 c0 8000 critical length L0 7500 7000 2 4 6 8 10 12 Axon Length (100 µm) 14 Figure 5.14. Relationship of the mean protein output c̄ and axonal length L, obtained by time averaging the solution to equation (5.34) for several values of τ. Function definitions and parameter values are as in Figure 5.13. The existence of a threshold protein output c0 could provide a mechanism for determining a critical length L0 . unless it is reduced by a considerable amount. Thus, the mean protein output of the system is relatively insensitive to the amplitude of the input signal and responds only to the frequency of the input signal, making the feed forward serial network a plausible means by which a neuron can decode the oscillating retrograde signal from the delayed feedback model. The monotonic relationship between c̄ and L suggests that the underlying intrinsic axonal length sensor could be based on a threshold protein value. That is, suppose that a given neuron is preprogrammed to grow until the mean protein output reaches some threshold value, c0 , see Figure 5.14. Based on the mean protein output, the neuron would be able to sense its critical length L0 and stop growing, for example. Analogous thresholding mechanisms have been investigated within the context of intracellular protein concentration gradients, which are used to determine spatial position within a cell so that, for example, cell division occurs at the appropriate time and location [116, 118]. Similarly, developmental morphogen gradients control patterns of gene expression so that each stage of cell differentiation occurs at the correct spatial location within an embryo. For biological effectiveness, these gradient-based mechanisms must be robust to intrinsic and extrinsic cellular noise [49, 116]. The issue of robustness to noise carries over to the 116 proposed axonal-length sensing mechanism, and can be analyzed along similar lines to protein concentration gradients. Therefore, we now investigate the impact of intrinsic noise in a gene network arising from finite copy numbers on the shape of the deterministic c̄ versus L curve. 5.5 Effects of Intrinsic Noise on Axonal Length Sensing In order to investigate the effects of intrinsic noise, we consider a slightly modified version of the network analyzed in section 5.5. Suppose that a gene promoter has two states: an inactive state Q and an active state Q∗ . In the active state, the gene produces the protein C at a rate of µ, and the protein subsequently decays at a rate λ. The promoter is activated in response to the pulsatile retrograde signal u I (t) and deactivates at a constant rate of β; see Figure 5.15: u I (t) µ λ −− Q) →C− → ∅. −* − Q∗ − (5.40) β Suppose there are N total gene promoters, each of which can exist in an active state or an inactive state. If N is sufficiently large, then the effects of intrinsic noise are negligible and one can represent the deterministic dynamics using kinetic equations. Let p(t) and x (t) denote, respectively, the fraction of active genes and the concentration of proteins (number of proteins per gene) at time t. Then dp = s(t)(1 − p) − βp, dt dx = µp(t) − λx (t), dt (5.41) where s(t) = u I (t) is the oscillatory retrograde signal coming from the delayed feedback model. Assume without loss of generality that p(0) = 0, so that the solution takes the form Z t Z z p(t) = s(z)exp β(z − t) + s(ξ )dξ dz. (5.42) 0 t We would like to calculate the time-averaged level of active genes in the large-time limit. In order to simplify our calculations, we proceed as in section 5.5 and take the oscillatory signal to consist of square pulses of width η and period T. Setting t = MT, positive integer M, we can break up the integral on the right-hand side of equation (5.42) into a sum of integrals evaluated over a single period: Z M −1 Z ( n +1) T p( MT ) = ∑ s(z)exp β(z − MT ) + n =0 = nT z MT s(ξ )dξ dz i exp(− βMT − Mη ) − 1 1 h exp ( β + 1)η − 1 , β+1 1 − exp( βT + η ) (5.43) (5.44) 117 Transcription factor activity Oscillatory signal T Promoter switching Protein output uI(t) Q nucleus time t uI(t) cytoplasm λ β Q* μ Figure 5.15. A gene promoter driven by the oscillatory retrograde signal u I (t). Adapted and redrawn from [117]. The second line comes from evaluating the various integrals and summing the resulting geometric series. Taking the limit M → ∞ shows that p( MT ) → Γ with Γ≡ i 1 h 1 . 1 − exp ( β + 1)η β+1 1 − exp( βT + η ) (5.45) For t ∈ [0, η ], we have dp0 1 1 − e−( β+1)t + Γe−( β+1)t , = 1 − (1 + β ) p0 ⇒ p0 ( t ) = dt 1+β (5.46) whereas for t ∈ (η, T ], dp1 = − βp1 ⇒ p1 (t) = p0 (η )e− β(t−η ) . dt (5.47) We have imposed continuity of the solution at t = η. Finally, p̄ is obtained by averaging the resulting periodic function over [0, T ]. p̄ = = 1h T Z η 0 p0 dt + Z T p1 dt η i i 1 1 1 − e β ( T − η ) η ( β +1) − 1 ) , + ( e β+1 T Tβ( β + 1) 1 − e βT +η hη (5.48) The formula for p̄ given above can be intuited in the following way. The fraction η/T corresponds to the fraction of time that s(t) is "on". The latter term in the bracketed sum is a correction for the alterations in the time-scale of the gene promoter reaction to s(t). When s(t) = 1, the time-scale of the gene promoter response is given by ( β + 1)−1 , whereas when s(t) = 0, the time-scale is given by β−1 . Finally, the time-averaged protein output is c̄ = µ p̄. λ (5.49) As in the simpler gene network of section 5.5, the time-averaged protein output is a monotonically decreasing function of T. 118 Now suppose that N is sufficiently small so that fluctuations due to low copy numbers cannot be ignored. In order to calculate the size of fluctuations, we have to consider the chemical master equation of the reaction scheme (5.40). Let n1 denote the total number of activated genes and let n2 denote the number of proteins that are present. Let P ≡ P(n1 , n2 , t) denote the probability that at a given time t, there are n1 active genes and n2 proteins available. The master equation is then given by dP =s(t)( N − n1 + 1) P(n1 − 1, n2 , t) + β(n1 + 1) P(n1 + 1, n2 , t) dt + µn1 P(n1 , n2 − 1, t) + λ(n2 + 1) P(n1 , n2 + 1, t) − (s(t)( N − n1 ) + βn1 + µn1 + λn2 ) P(n1 , n2 , t). (5.50) The first two terms correspond to the activation or the deactivation of a gene that results in having n1 active genes and n2 proteins. The second two terms correspond to the production or the degradation of a protein that results in having n1 active genes and n2 proteins. The last terms correspond to the ways that the system can leave the state of having n1 active genes and n2 proteins. It is difficult to solve the master equation explicitly, so we carry out a system size expansion with respect to N. That is, set ni = Nxi and rewrite equation (5.50) as dP 1 1 = N [s(t)(1 − x1 + ) P(n1 − 1, n2 , t) + β( x1 + ) P(n1 + 1, n2 , t) dt N N 1 + λx1 P(n1 , n2 − 1, t) + γ( x2 + ) P(n1 , n2 + 1, t) N − (s(t)(1 − x1 ) + βx1 + λx1 + γx2 ) P(n1 , n2 , t)]. (5.51) (5.52) (5.53) The master equation is now just a sum of terms of the form f (n/N ) P(n, t), where n ≡ (n1 , n2 ) and f is the corresponding propensity function. Performing the change of variables f (n/N ) P(n, t) = f (x) p(x, t), where x ≡ ( x1 , x2 ), and Taylor expanding in powers of N −1 to second-order leads to the Fokker-Planck equation dp ∂ ∂ =− s(t)(1 − x1 ) − βx1 − (µx1 − λx2 ) p dt ∂x1 ∂x2 1 ∂2 ∂2 + ( s ( t )( 1 − x ) + βx ) p + ( µx + λx ) . 2 1 1 1 2N ∂x12 ∂x22 (5.54) (5.55) In the case of a constant input s(t) = α, the deterministic kinetic equations (5.41) have the unique fixed point solution x1∗ = α , α+β x2 ∗ = µ ∗ x . λ 1 (5.56) 119 In this case, neglecting transients, the Fokker-Planck equation describes a stochastic process characterized by Gaussian fluctuation about the fixed point ( x1∗ , x2∗ ). It is then relatively straightforward to calculate the stationary variance ∆c of the protein output, given the mean hci = x1∗ : αµ i µ2 αβ 1h + γ (α + β + λ)(α + β)2 α+β h i µ α µβ = 1+ λα+β (α + β + λ)(α + β) h i µβ = hci 1 + . (α + β + λ)(α + β) ∆c = (5.57) (5.58) (5.59) The expression for the variance in the case of constant input consists of an intrinsic Poissonian term due to random protein production and an extrinsic term due to fluctuations in the gene promoters themselves. The calculation of the variance in the case of an oscillatory input s(t) is considerably more involved, even when it takes the form of square pulses. However, stochastic simulations show that the protein variance in response to an oscillatory signal is less than the protein variance in response to a constant input, assuming that time-averaged means are the same [117]. Let αeff be the effective constant input for which the time-averaged and noise-averaged protein output hc̄i can be written as hc̄i = µ αeff . λ αeff + β (5.60) It follows that for an oscillatory input h ∆c̄ ≤ hc̄i 1 + i µβ , (αeff + β + λ)(αeff + β) (5.61) where ∆c̄( L) = h[c̄( L) − hc̄( L)i]2 i. 5.5.1 (5.62) Errors in Axonal Length Sensing The relationship between c̄ and L in the presence of intrinsic noise is shown in Figure 5.16. The general inverse relationship is still prevalent in this situation, but fluctuates due to the stochasticity in the gene switching. By analogy with the effects of intrinsic noise in protein concentration gradients [49], the presence of noise in the protein output leads to 120 9800 0.4 (a) 0.35 9700 Mean Protein Output 0.3 ∆ L (100 µm) 9600 0.25 9500 0.2 0.15 9400 0.1 9300 9200 (b) 0.05 2 4 6 Axon Length (100 µm) Relative Error (∆ L/LT) 0.045 0.04 8 0 10 2 4 6 8 Threshold Length (100 µm) 10 (c) 0.035 0.03 0.025 0.02 0.015 0.01 2 4 6 8 Axonal Length (100 µm) 10 Figure 5.16. Results of simulation of the chemical master equation (5.50). (a) Plot of mean protein output c̄ versus axonal length L. (b) Plot of uncertainty in axonal length ∆L versus threshold axonal lengths L T . (c) Relative error (∆L/L T ) versus axonal length. Parameter values used to generate retrograde signal u I (t) are the same as in Figure 5.13. Other parameter values are β = 1, µ = .1, λ = 0.01, and N = 1000. an uncertainty ∆L in the critical axonal length L at which the threshold is crossed. This is illustrated schematically in Figure 5.17. The uncertainty ∆L can be estimated using q 0 ∆L|hc̄ ( L)i| = ∆c̄( L), (5.63) If we ignore the correction factor in equation (5.61) and approximate the stochastic process by a Poisson process, then ∆c̄ ≈ c̄ and ∆L ∼ p hc̄( L)i . |hc̄0 ( L)i| (5.64) As a further approximation, suppose that hc̄( L)i ∼ 1/T, where T is the period of oscillations produced by an axon of length L, so that |hc̄0 ( L)i| ∼ 1/( L0 ( T ) T 2 ) (using the fact that _ mean protein output c(L) 121 _ <c> _ Δc cT ΔL axonal length L Figure 5.17. Results of simulation of the chemical master equation (5.50). (a) Plot of mean protein output c̄ versus axonal length L. (b) Plot of uncertainty in axonal length ∆L versus threshold axonal lengths L T . (c) Relative error (∆L/L T ) versus axonal length. Parameter values used to generate retrograde signal u I (t) are the same as in Figure 5.13. Other parameter values are β = 1, µ = .1, λ = 0.01, and N = 1000. L increases monotonically with T). It follows that ∆L T 3/2 L0 ( T ) ∼ , L L( T ) (5.65) Assuming that the length increases at least linearly with T, we see that the relative error grows with the oscillation period T and, hence, the axonal length L. Although this is a crude estimate, we find that the same qualitative behavior is observed in numerical simulations of the full stochastic model. This is shown in Figure 5.16c, where we plot the relative error ∆L/L T versus axonal length. Our analysis suggests that the frequencyencoded protein threshold mechanism could break down for long axons. An analogous result was shown to hold in [56], where the robustness of the encoding of axonal length in the frequency of a pulsatile signal was investigated. There we found that the encoding of axonal length into frequency became less reliable at long axon lengths due to accumulation of white noise signified by a high coefficient of variation in the frequency of the retrograde signal. In this work, the retrograde signal is deterministic, and the error in protein output 122 is accounted for strictly by the random variations in the activities of independent gene promoters. Hence, the error in length sensing could be more devastating in real life situations, since noise would impact both the encoding and the decoding processes. Thus, wherever the sources of noise may be, their impact on this frequency-dependent mechanism is clear: large neurons would have a more difficult time sensing their own length when compared with smaller neurons. 5.6 Discussion In this chapter, we developed a mathematical model for axonal length sensing based on a biophysical mechanism recently proposed by Rishal et al. [102]. We showed that the underlying dynamical mechanism involves delayed negative feedback due to the finite propagation speeds of molecular motors. This can be incorporated into the kinetic equations for retrograde chemical signaling using a discrete delay τ = L/v or by convolving the chemical signals with the Green's function of an advection-diffusion equation for motor transport. Both versions of the model support chemical oscillations that emerge via a Hopf bifurcation, resulting in a frequency that is inversely related to the axonal length. Furthermore, the advection-diffusion version of the model suggests that knockdown of either kinesin or dynein motors results in a longer axon (see Figure 5.9). These results are consistent with the experimental and computational studies carried out in [102]. The advantage of our mathematical model is that it provides a compact dynamical framework for understanding the origin of the oscillations and for exploring how the length-sensing mechanism depends on various biophysical parameters. One prediction of our model is that the effective diffusivity D of motor transport only has a weak affect on the retrograde signal frequency's dependence on axonal length. That is, increasing D increases the critical length Lc for the onset of oscillations in the retrograde signal, but once oscillations have formed, the frequency is approximately D-independent. A second prediction is that fluctuations in the transport of molecular motors results in a reduction in the reliability of the frequency-encoding mechanism as the length increases (see Figure 5.11). This could have significant implications for the viability of such a mechanism within the context of axonal injury, where accurate information regarding the location of the injury is needed in order to target regeneration at the correct location [55]. 123 We also showed that the mean protein output of a gene network responding to the pulsatile retrograde signal varies inversely with axonal length. Specifically, we derived approximate analytical results which make explicit the inverse relationship, and introduced the notion that frequency decoding could be done based on a protein threshold mechanism. We then investigated the reliability of such a mechanism subject to intrinsic noise stemming from finite copy numbers within a gene network by numerically simulating a chemical master equation describing the random switching of genes and production of protein. The results of these simulations suggest that the accuracy in the information the neuron receives regarding axonal length declines as axonal length itself grows. This result could have serious implications for the utility of this mechanism in the context of axonal injury, where accurate information regarding the locality of an affliction is necessary for a neuron to set a regenerative process in motion. There are a number of possible extensions of our model of axonal length control. First, it would be interesting to derive an effective stochastic advection-diffusion equation with multiplicative noise, starting from a full stochastic model of stop-and-go motor transport. Such a model would need to keep track of the total number of kinesin and dynein motors. In our simplified model, we assumed that there was a sufficient supply of kinesin motors at the proximal end and dynein motors at the distal end. The model was self-regulating due to the feedback signals. Secondly, equation (5.39) indicates that mean protein output is sensitive to the frequency of the incoming pulse signal, but that it also sensitive to the fraction of time for which the incoming signal is at its peak value. Hence, the feed forward serial network does not generate a frequency filter in the strictest sense. We would be interested in seeing the result of feeding the retrograde signal into a network that allows for more acute frequency sensitivity. Yet another extension of our work here would be to inject kinesin motors into the axon based on a doubly stochastic Poisson process as in Chapter 4 for flagellar length control. This would incorporate another source of noise that would impact the frequency-encoding of length information of axons. Furthermore, we have represented an axon as a 1D domain when in reality axons are highly branched in structure. Considering a tree-type domain for this process would be more appropriate. Finally, it would be interesting to build off of our work in Chapter 3 and incorporate exclusion effects in our axon-length sensing model. CHAPTER 6 FUTURE DIRECTIONS In this dissertation, the mathematical models used to describe complex biological processes involved the invocation of particular approximations that facilitated mathematical analysis. For example, as discussed in Chapter 1, we assume that cargo elements are carried by a single motor whose dynamics may be described by a three-state model, and in Chapter 2 we assume that the distribution of microtubules in two- and three-dimensional cells is symmetric. The descriptive limitations of our models have therefore motivated several problems to be worked on in the future. 6.1 Microscopic Descriptions of Motor Transport In Chapter 1, we discussed a mathematical formulation of a tug-of-war mechanism for bidirectional motion of cargo within a cell; see equations (1.8) and (1.9). To facilitate mathematical analysis, we again have to invoke an adiabatic approximation on equation (1.6) with the matrices and vectors formatted to reflect the tug-of-war mechanism. The resulting Fokker-Planck equation's drift and diffusion terms would be significantly more complicated, but, provided that we model a large population of motor-cargo complexes evolving according to this equation, we may use it as the dictating partial differential equation for bulk motor dynamics. We may then couple the FP equation with processes such as vesicular transport and size homeostasis. The adiabatic reduction process in itself is more complicated than that of the three-state model; solving equation (1.19), for example, involves the use of the singular value decomposition (SVD) [18]. The drift and diffusive terms are functions of various biophysical parameters of the tug-of-war model. These include the stall force Fs , the detachment force Fd , the maximum forward and backward velocities, v f , vb , and binding/unbinding rates. Hence, we may be able to determine on a microscopic level the critical parameters underlying a cell's ability to achieve uniformity in vesicle distribution across its body. 125 Furthermore, it is unclear whether motors participating in tug-of-war type dynamics are disparate but located on the same microtubule or the same but situated on adjacent microtubules of opposite polarity. This idea of navigating between several microtubules nearby should also come into play when thinking about exclusion effects. In our work, we have assumed that motors carrying cargo may not be passed by other motors due to exclusion effects; however, it is entirely possible that motors may hop to an adjacent microtubule and move past other motor-cargo complexes. Furthermore, as discussed in Chapter 1, microtubules are large structures, and there is plenty of room for motors to move past one another. Work in Chapter 3 was motivated by experimental observations of molecular motor jamming [86]. Hence, future work should entail models that account for networks of microtubules, potentially of differing polarity, along which motors can move. This type of modeling would be relevant for our work on axonal length sensing as well. Another central aspect of active transport that needs to be addressed is the individual chemical reactions that underly kinesin and dynein dynamics. The work in this dissertation largely ignores the microscopic reactions that involve the use of ATP to propel kinesin and dynein motors forward. Furthermore, we have assumed that motor interaction with vesicles is governed by first-order kinetics. Recent work has shown that these kinetics may be driven by Michaelis-Menten-type interactions [23]; the equations in Chapter 2 and Chapter 3 are more complicated and may lead to more interesting mathematics in looking at the steady-state distribution of vesicles. Accounting for these reactions and modeling cargo dynamics with a tug-of-war mechanism may lead to more interesting biophysical insights regarding the crucial underlying parameters that determine whether or not uniformity in vesicle density may be achieved. 6.2 Motor Routing by Microtubule Network Modification There is now growing evidence that the routing of molecular motors is achieved by modifying the shape of networks of microtubules [127]. This is a new wrinkle in the question of how cells route motors to delivery specific cargo to specific locales in a cell. This involves production of protein involved in a signal pathway microtubule length, (-) end spacing, and coverage. It would be interesting to look at the phenomenon of uniformity in vesicle distribution from the perspective of microtubule networks and the cell's nucleus 126 rather than from the perspective of motor-cargo interactions. The nature of the modeling would seemingly have a multiscale flavor to it. To account for protein production, one must study the state of the particular gene network producing the said protein. Then the nature of the pathway involving the protein altering microtubule network structure must be examined; furthermore, the mechanism of gene regulation via feedback must be considered. It would be interesting to see, from this perspective, what biophysical parameters facilitate, for example, uniform vesicle distribution across a given cell's interior. 6.3 Neuron Polarization In section 4.4, we briefly mention that an extension of the work in Chapter 4 is to consider cases where two flagella are competing for resources from some finite pool to outgrow one another. A similar idea may be used to investigate the genesis of axons in nascent neurons. Neurons are highly polarized cells, typically having a single long, thin axon and several short, thick dendrites. In a canonical neuron, the axon transmits information to the neuron's target and dendrites obtain and process information. The question of how this polarization is achieved, however, has been the subject of intense examination over the past 40 years. It is clear that neuron polarization occurs during early stages of cell differentiation, before migration [20]. Neuron polarization appears to undergo a 5-step procedure [26]. In stage 1 of polarization, neurons extend a motile lamellipodia around the cell body. In stage 2, the lamellipodia clusters at particular, relatively symmetrical sites until neurites form. These neurites are highly dynamic, growing and retracting in a stochastic fashion until one of them undergoes a sudden and sustained growth; this singular neurite develops into the axon. The formation of the axon characterizes stage 3. In stage 4, the remaining neurites develop into dendrites. In stage 5, synaptic specializations and contacts are established [20]. The fact that neuronal polarity is apparent by stage 3 suggests that polarity is established in the stage-2-stage-3 transition. Researchers have been searching for a molecular basis for the establishment of this polarity, and several theoretical models for how these molecules establish polarity have been proposed (see, for example, [5, 33, 40, 44, 84, 115]). In Ref. [115], a mathematical model for spontaneous neuronal symmetry breaking was 127 developed based on experimentally observed dynamics of a molecule known as shootin-1 during the stage-2-stage-3 transition. Stage-2 neurons are relatively symmetric and one characteristic of this phase is the stochastic accumulation and dissipation of shootin-1 at the growth cones of each of the neurites, which all have approximately the same length. Eventually a significant proportion of the shootin-1 in a neuron accumulates in one of the neurite growth cones and causes an immediate outgrowth of that neurite; this becomes the axon. The neuronal outgrowth is caused by the accumulation of shootin-1 in the growth cone inducing a traction force. The traction force is opposed by a neurite length-dependent tension. Once one of the neurites has been selected to be the axon, it accumulates more and more of the available shootin-1. Shootin-1 is transported from the cell body to the growth cone in packed boluses via anterograde active transport; retrograde transport is achieved by means of passive diffusion. Hence, lengthier neurites keep shootin-1 in their interior longer than shorter ones. The result is a positive feedback loop: (i) accumulation of shootin-1 induces neurite outgrowth and (ii) neurite outgrowth promotes shootin-1 to stay in its interior longer. We can capture the "big picture" dynamics of the above-described mechanism by considering the concentration of shootin-1 in the growth cone of a given neurite and the length of that neurite. The dynamics of the shootin-1 concentration are AD dC =− (C − C0 ) + w(t), dt VL (6.1) where C (t) is the concentration of shootin-1 in the growth cone, C0 is the concentration of shootin-1 in the cell body, A is the cross-sectional area of the neurite, D is the diffusion coefficient, V is the volume of the growth cone, L is the length of the neurite, and w(t) is a function describing the arrival of shootin-1 to the growth cone. It should reflect that the arrival times are stochastic. An example is w(t) = wavg ∑ g(t − t j ), (6.2) j where wavg is the average concentration of shootin-1 in each bolus, t j is a stochastic arrival time, and g(t) is a Gaussian function. The dynamics of the length are determined by the traction force induced by shootin-1 accumulation and neurite length-dependent tension. Hence, dL = δk on M − δk o f f exp − ( F (C ) − T ( L)) , dt (6.3) 128 where M is the concentration of the material used to construct the neurite shaft (typically tubulin), k on,o f f are the rates at which the units of the construction material bind and unbind from the neurite shaft, δ is the length of each of these units, F (C ) is the traction force, and T ( L) is the tension. In Ref. [115], a model for neuron polarization were constructed from equations similar to equations (6.1) and (6.3). Parameter values and equations for the forces were obtained by fitting to actual experimental data. Numerical simulations showed that their model was sufficient for neuron symmetry-breaking. Furthermore, their model showed that such a mechanism was sufficient for axon regrowth after transection, which is observed experimentally. Even more incredibly, their model showed cases where neurons developed multiple axons, which is observed in some experiments as well. It would be interesting to develop these types of models from the perspective of molecular motors. For example, instead of including a source of shootin-1 with a function w(t) as in equation (6.1), we could model anterograde motion of shootin-1 down a neurite with an advection-diffusion equation and encode the shootin-1 concentration, C (t), into the boundary conditions. Similarly, we could model retrograde motion with a diffusion equation. We could then include additive noise in our model as we did in Chapter 5 for our axonal length-sensing model or bring in ideas from Chapter 4 and allow the injection times of shootin-1 into the neurite to be determined by a doubly stochastic Poisson process. 6.4 Closing Remarks As all models are limited, ours fail to address the aforementioned aspects of important processes intertwined with active intracellular transport. But we have helped open the door to mathematical modeling of size homeostasis and vesicular delivery. The deep insights of new research suggests that the full picture of how active transport works may only be understood upon constructing a multiscale model of active transport and its regulation. These prospects are very exciting. APPENDIX GREEN'S FUNCTION FOR ADVECTION DIFFUSION EQUATION In this appendix, we derive the explicit formulae (5.18a) and (5.18b) in Chapter 5 for the Green's functions G− (ξ, σ; x, t) of the advection-diffusion model. For concreteness, we will focus on the Green's function G+ , since the derivation of G− is very similar. The Green's function G+ is defined according to the equation − ∂G+ ∂2 G+ ∂G+ = δ ( x − ξ ) δ ( t − σ ), −v −D ∂σ ∂ξ ∂ξ 2 where δ( x ) is the Dirac-Delta function, supplemented by the boundary conditions vG+ (ξ, σ; x, t) + D ∂G+ (ξ, σ; x, t) = 0, ∂ξ ξ = 0, L (A.1) Introduce the linear differential operator L acting on functions u = u(ξ, σ), Lu ≡ ∂u ∂u ∂2 u +v −D 2, ∂σ ∂ξ ∂ξ and define the inner product h f , gi ≡ Z LZ ∞ 0 f g dtdx. 0 We can then rewrite the equation for G in the more compact form L† G+ (ξ, σ; x, t) = δ( x − ξ )δ(t − σ), (A.2) where L† is the adjoint operator L† u ≡ − ∂u ∂u ∂2 u −v −D 2. ∂σ ∂ξ ∂ξ Using integration by parts and the boundary conditions on G and the kinesin concentration c+ , one can show that h G+ , Lc+ i = hc+ , L† G+ i − − Z ∞ 0 Z L 0 G+ (ξ, 0; x, t)c+ (ξ, 0)dξ G+ (0, σ; x, t)v ( IE − w I f [u I (σ)]) dσ, 130 which reduces to equation (5.16), since h G+ , Lc+ i = 0, hc+ , L† G+ i = c+ ( x, t) and G+ (0, σ; x, t) = 0 for σ > t (causality). It remains to solve the adjoint problem given by equation (A.2). Introduce the change of variables σ̄ = t − σ and set G+ (ξ, σ; x, t) = G (ξ; x, σ̄). We then have ∂G ∂2 G ∂G −v − D 2 = δ(ξ − x )δ(σ̄) ∂σ̄ ∂ξ ∂ξ (A.3) with G (ξ; x, σ̄) = 0 for σ̄ < 0. Applying the Laplace transform to equation (A.3) with e(ξ; x, s) = G we obtain Z ∞ 0 e−st G (ξ; x, t)dt, e − vG e0 − D G e00 = δ(ξ − x ) sG where 0 denotes differentiation with respect to ξ for fixed s, x. It is convenient to eliminate e(ξ; x, s) = g(ξ; x, s)φ(ξ ) for an appropriately chosen the first derivative term by setting G e gives function φ. Substituting into the equation for G − Dφg00 − (vφ + 2Dφ0 ) g0 + (sφ − vφ0 − Dφ00 ) g = δ(ξ − x ) v v Hence, imposing φ0 = − 2D φ ⇒ φ = e(− 2D ξ ) , we see that g satisfies the self-adjoint (Sturm- Liouiville) equation g00 + − v2 + 4Ds 1 v g = − e( 2D x) δ(ξ − x ) 2 4D D (A.4) We can then solve this equation in terms of the eigenfunctions ζ n and eigenvalues λn of the second-order equation ζ n00 + λn ζ n = 0, supplemented by the homogeneous boundary conditions v dζ n − ζ n (0) + D (0) = 0 2 dξ (A.5) v dζ n ζ n ( L) + D ( L) = 0. 2 dξ (A.6) The latter follow from the boundary condition for G+ . Using the fact that the eigenfunctions form a complete orthonormal set, we have the expansions ∞ δ(ξ − x ) = ∑ ζ n ( x )ζ n (ξ ) n =1 131 and ∞ g(ξ; x, s) = ∑ an (x, s)ζ n (ξ ) (A.7) n =1 Substituting these into (A.4), we can solve for an ( x, s) to obtain ∞ g(ξ; x, s) = ∑ n =1 1 s + Dλn + v v2 4D e( 2D x) ζ n ( x )ζ n (ξ ) (A.8) e inverting the Laplace transform, and Substituting g and φ back into the formula for G, reverting back to original time coordinates, we finally obtain ∞ G+ (ξ, σ; x, t) = ∑ e(−(Dλ + n v2 4D )( t − σ )) v e(− 2D (ξ − x)) ζ n (ξ )ζ n ( x ) n =1 Using standard methods to solve boundary value problems, we find that √ p p 2D λ ζ n = sin λn x + cos ( λn x ), v where λn solves 4Dv p We thus obtain equation (5.18a). λn cot ( p λn L) = 4D2 λn − v2 . (A.9) REFERENCES [1] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, Molecular Biology of the Cell, Garland, 5th ed., 2008. [2] T. Alberts, K. Khanin, and J. Quastel, The continuum directed random polymer, Journal of Statistical Physics, 154 (2014). [3] C. A. Albus, I. Rishal, and M. Fainzilber, Cell length sensing for neuronal growth control., Trends Cell Biol., 23 (2013), pp. 305-310. [4] C. Appert-Rolland, J. Cividini, and H. J. Hilhorst, Intersection of two TASEP traffic lanes with frozen shuffle update, Journal of Statistical Mechanics: Theory and Experiment, 10 (2011). [5] N. Arimura and K. Kaibuchi, Neuronal polarity: from extracellular signals to intracellular mechanisms, Nat. Rev. Neurosci., 8 (2007). [6] M. Basu and P. K. Mohanty, Asymmetric simple exclusion process on a Cayley tree., Journal of Statistical Mechanics: Theory and Experiment, 10 (2010). [7] J. Best, Doubly stochastic processes: an approach for understanding central nervous system activity., Proc. 3rd Int. Conf. Appl. Math. Modeling, (2009), pp. 155-158. [8] R. A. Blythe and M. R. Evans, Nonequilibrium steady steady states of matrix-product form: a solver's guide, J. Phys. A., 40 (2007), pp. R333-R441. [9] A. Borodin, I. Corwin, and T. Sasamoto, From duality to determinants for q-TASEP and ASEP, The Annals of Probability, 42 (2014), pp. 2314-2382. [10] P. R. Bouzas, M. J. Valderrama, and A. M. Aguilera, On the characterisitc functional of a doubly stochastic Poisson process: Application to a narrow-band process., Appl. Math. Modelling, 30 (2006), pp. 1021-1032. [11] P. C. Bressloff, A stochastic model of intraflagellar transport, Phys. Rev. E, 73 (2006). [12] P. C. Bressloff, Stochastic Processes in Cell Biology, Springer, 2014. [13] P. C. Bressloff, Aggregation-fragmentation model of vesicular transport in neurons, J. Phys. A., 49 (2016). [14] P. C. Bressloff and B. R. Karamched, A frequency-dependent decoding mechanism for axonal length sensing, Front. Cell. Neurosci., 9 (2015). [15] P. C. Bressloff and B. R. Karamched, Doubly stochastic Poisson model of flagellar length control, Submitted., (2016). [16] P. C. Bressloff and B. R. Karamched, Model of reversible vesicular transport with exclusion, J. Phys. A., 49 (2016). 133 [17] P. C. Bressloff and E. Levien, Synaptic democracy and active intracellular transport in axons, Phys. Rev. Lett., 114 (2015). [18] P. C. Bressloff and J. Newby, Stochastic models of intracellular transport, Rev. Mod. Phys., 85 (2013), pp. 135-196. [19] P. C. Bressloff and J. M. Newby, Quasi-steady state analysis of motor-driven transport on a two-dimensional microtubular network., Phys. Rev. E, 83 (2011). [20] A. Caceres, B. Ye, and C. G. Dotti, Neuronal polarity: demarcation, growth and commitment, Current Opinion in Cell Biology, 24 (2012), pp. 547-553. [21] L. Cai, C. K. Dalal, and M. B. Elowitz, Frequency-modulated nuclear localization bursts coordinate gene regulation, Nature, 45 (2008), pp. 485-491. [22] T. Chou, K. Mallick, and R. K. Zia, Non-equilibrium statistical mechanics: from a paradigmatic model to biological transport, Rep. Prog. Phys., 74 (2011). [23] L. Ciandrini, M. C. Romano, and A. Parmeggiani, Stepping and crowding of molecular motors: Statistical kinetics from an exclusion process perspective, Biophys. J., 107 (2014), pp. 1176-1184. [24] R. Corless, G. Gonnet, D. Hare, D. Jeffrey, and D. Knuth, On the Lambert W function, Adv. Comput. Math., 5 (1996). [25] D. R. Cox, Some statistical methods connected with series of events., Journal of Royal Statistical Society B, 17 (1955), pp. 129-164. [26] A. M. Craig and G. A. Banker, Neuronal polarity., Annu. Rev. Neurosci., 17 (1994), pp. 267-319. [27] J. A. Deane, D. G. Cole, and J. L. Rosenbaum, Localization of intraflagellar transport protein IFT 52 identifies basal body transitional fibers as the docking site for IFT particles., Curr. Biol., 11 (2001), pp. 1586-1590. [28] B. Derrida, An exactly solvable non-equilibrium system: the asymmetric simple exclusion process., Phys. Rep., 301 (1998). [29] B. Derrida, E. Domany, and D. Mukamel, An exact solution of a one-dimensional asymmetric exclusion model with open boundaries., J. Stat. Phys., 69 (1992). [30] B. Derrida, M. R. Evans, V. Hakim, and V. Pasquier, Exact solution of a 1D asymmetric exclusion model using a matrix formulation., J. Phys. A., 26 (1993). [31] D. W. Erickson, G. Pruessner, B. Schmittmann, and R. K. P. Zia, Spurious phase in a model for traffic on a bridge, J. Phys. A., 38 (2005). [32] M. R. Evans, R. Juhasz, and L. Santen, Shock formation in an exclusion process with creation and annihilation, Phys. Rev. E, 68 (2003). [33] M. Fivaz, S. Bandara, T. Inoue, and T. Meyer, Robust neuronal symmetry breaking by ras-triggered local positive feedback, Curr. Biol., 18 (2008), pp. 44-50. [34] C. W. Gardiner, Handbook of Stochastic Processes., Springer, Berlin., 2009. 134 [35] D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions., The Journal of Physical Chemistry, 81.25 (1977), pp. 2340-2361. [36] N. W. Goehring and A. A. Hyman, Organelle growth control through limiting pools of cytoplasmic components., Current Biology, 22 (2012), pp. R330-R339. [37] J. L. Goldberg, How do axons grow., Genes and Dev., 17 (2009), pp. 941-958. [38] K. Goslin and G. Banker, Experimental observations on the development of polarity by hippocampal neurons in culture, J. Cell Biol., 108 (1989), pp. 1507-1516. [39] B. S. Govindan, M. Gopalakrishnan, and D. Chowdhury, Length control of microtubules by depolymerizing motor proteins., Europhys. Lett., 83 (2008). [40] B. P. Graham and A. van Ooyen, Mathematical modelling and numerical simulation of the morphological development of neurons, BMC Neuroscience, 7(Suppl 1) (2006). [41] J. Grandell, Doubly Stochastic Process., Springer-Verlang., 1976. [42] S. P. Gross, Hither and yon: a review of bi-directional microtubule-based transport, Phys. Biol., 1 (2004), pp. R1-11. [43] P. W. Gunning, U. Ghoshdastider, S. Whitaker, D. Popp, and R. C. Robinson, The evolution of compositionally and functionally distinct actin filaments, Journal of Cell Science, 128 (2015), pp. 2009-2019. [44] J. J. J. Hjorth, J. van Pelt, H. D. Mansvelder, and A. van Ooyen, Competitive dynamics during resource-driven neurite outgrowth, PLoS ONE, 9 (2014). [45] F. J. Hoerndli, D. A. Maxfield, P. J. Brockie, J. E. Mellem, E. Jensen, R. Wang, D. M. Madsen, and A. V. Maricq, Kinesin-1 regulates synaptic strength by mediating the delivery, removal, and redistribution of AMPA receptors, Neuron, 80 (2013), pp. 1421- 1437. [46] W. Hong, A. Takshak, O. Osunbayo, A. Kunwar, and M. Vershinin, The effect of temperature on microtubule-based transport by cytoplasmic dynein and kinesin-1 motors., Biophys J., 111 (2016), pp. 1287-1294. [47] L. E. Hough, A. Schwabe, M. A. Glaser, J. R. McIntosh, and M. D. Betterton, Microtubule depolymerization by kinesin-8 motor kip3p: a mathematical model., Biophys. J., 96 (2009), pp. 3050-3064. [48] J. Howard, Mechanics of Motor Proteins and the Cytoskeleton, Sinauer, 2001. [49] M. Howard, How to build a robust intracellular concentration gradient, Trends Cell Biol., 22 (2012), pp. 311-317. [50] H. S. Huan and M. D. Betterton, Biophysics of filament length regulation by molecular motors., Phys. Biol., 10 (2013). [51] F. Huber, J. Schnauss, S. Ronicke, P. Rauch, K. Muller, C. Futterer, and J. Kas, Emergent complexity of the cytoskeleton: from single filaments to tissue, Advances in Physics, 62 (2013), pp. 1-112. 135 [52] T. Imamura and T. Sasamoto, Stationary correlations for the 1D KPZ equation, Journal of Statistical Physics, 150 (2013), pp. 908-939. [53] D. Johann, C. Erlenkamper, and K. Kruse, Length regulation of active biopolymers by molecular motors, Phys. Rev. Lett., 108 (2012). [54] M. Kac, On the distribution of certain Wiener functionals., Trans. Am. Math. Soc., 65 (1949), pp. 1-13. [55] N. Kam, Y. Pilpel, and M. Fainzilber, Can molecular motors drive distance measurements in injured neurons?, PLoS Comput. Biol., 5 (2009), p. e1000477. [56] B. R. Karamched and P. C. Bressloff, Delayed feedback model of axonal length sensing, Biophys. J., 108 (2015), pp. 2408-2419. [57] B. R. Karamched and P. C. Bressloff, Effects of cell geometry on reversible vesicular delivery, J. Phys. A., In press. (2016). [58] I. Katsura, Determination of bacteriophage λ tail length by a protein ruler., Nature, 327 (1987), pp. 73-75. [59] J. P. Keener, How salmonella typhimurium measures the length of flagellar filaments., Bull. Math. Biol., 68 (2006), pp. 1761-1778. [60] J. P. Keener and J. Sneyd, Mathematical Physiology, Springer, 1998. [61] M. J. Kennedy and M. D. Ehlers, Organelles and trafficking machinery for postsynaptic plasticity., Ann. Rev. Neurosci., 29 (2006). [62] K. Keren, Z. Pincus, G. M. Allen, E. L. Barnhart, G. Marriot, A. Mogilner, and J. A. Theriot, Mechanism of shape determination in motile cells, Nature, 453 (2008), pp. 475-480. [63] S. Klumpp and R. Lipowsky, Cooperative cargo transport by several molecular motors, PNAS, 102 (2005), pp. 17284-17289. [64] M. P. Koonce and M. Samso., Of rings and levers: the dynein motor comes of age., Trends Cell Biol., 14 (2004), pp. 612-619. [65] D. C. Krakauer, K. M. Page, and S. Sealfon, Module dynamics of the GnRH transduction network, J. Theor. Biol., 218 (2002), pp. 457-470. [66] J. Krug, Boundary-induced phase transitions in driven diffusive systems., Phys. Rev. Lett., 67 (1991), pp. 1182-1185. [67] M. Krumin and S. Shoham, Generation of spike trains with controlled auto- and crosscorrelation functions, Neural Computation, 21 (2009), pp. 1642-1664. [68] C. Kural, H. Kim, S. Syed, G. Goshima, V. I. Gelfand, and P. R. Selvin, Kinesin and dynein move a peroxisome in vivo: a tug-of-war or coordinated movement?, Science, 308 (2005), pp. 1469-1472. [69] T. Lagache and D. Holcman, Effective motion of a virus trafficking inside a biological cell., SIAM J. Appl. Math., 68 (2008), pp. 1146-1167. 136 [70] F. Lallemend and P. Ernfors, Molecular interactions underlying the specification of sensory neurons, Trends Neurosci., 35 (2012), pp. 373-381. [71] D. Lando, On Cox processes and credit risky securities., Review of Derivatives Research., 2 (1998), pp. 99-120. [72] S. D. Lawley, M. Tuft, and H. A. Brooks, Coarse-graining intermittent intracellular transport: Two- and three- dimensional models., Phys. Rev. E, 92 (2015). [73] P. A. Lefebvre and J. L. Rosenbaum, Regulation of the synthesis and assembly of ciliary and flagellar proteins during regeneration., Ann. Rev. Cell Dev. Biol., 2 (1986), pp. 517- 546. [74] M. Lemieux, S. Labrecque, C. Tardif, E. Labrie-Dion, E. Lebel, and P. D. Koninck, Translocation of CaMKII to dendritic microtubules supports the plasticity of local synapses., J. Cell Biol., 198 (2012), pp. 1055-1073. [75] P. A. Lewis and G. S. Shedler, Simulation of nonhomogeneous Poisson processes by thinning., Naval Research Logistics Quarterly., 26.3 (1979), pp. 403-413. [76] W. B. Ludington, H. Ishikawa, Y. V. Serebrenik, A. Ritter, R. A. HernandezLopez, J. Gunzenhauser, E. Kannegard, and W. F. Marshall, A systematic comparison of mathematical models for inherent measurement of ciliary length: how a cell can measure length and volume., Biophys. J., 108 (2015), pp. 1361-1379. [77] W. B. Ludington, K. A. Wemmer, K. F. Lechtreck, G. B. Witman, and W. F. Marshall, Avalanche-like behavior in ciliary import., Proc. Natl. Acad. Sci. USA, 110 (2013), pp. 3925-3930. [78] C. I. Maeder, A. San-Miguel, E. Y. Wu, H. Lu, and K. Shen, In vivo neuron-wide analysis of synaptic vesicle precursor trafficking, Traffic, 15 (2014), pp. 273-291. [79] S. N. Majumdar, Brownian functionals in physics and computer science., Curr. Sci., 89 (2005), pp. 2976-2092. [80] A. F. Maree, D. Jilkine, A. Dawes, V. A. Grieneisen, and L. Edelstein-Keshet, Polarization and movement of keratocytes: A multiscale modeling approach, Bull. Math. Biol., 68 (2006), pp. 1169-1211. [81] W. F. Marshall, H. Qin, M. R. Brenni, and J. L. Rosenbaum, Flagelar length control system: testing a simple model based on intraflagellar transport and turnover., Mol. Cell Biol., 16 (2005), pp. 270-278. [82] W. F. Marshall and J. L. Rosenbaum, Intraflagellar transport balances continuous turnover of outer doublet microtubules: implications for flagellar length control, J. Cell Biol., 155 (2001), pp. 405-414. [83] W. J. Marshall, Cellular length control, Ann. Rev. Cell Dev. Biol., 20 (2004), pp. 677- 693. [84] S. A. Menchon, A. Gartner, P. Roman, and G. Dotti, Neuronal (bi)polarity as a self-organized process enhanced by growing membrane, PLoS ONE, 6 (2011). 137 [85] P. Mottishaw, B. Waclaw, and M. R. Evans, An exclusion process on a tree with constant aggregate hopping rate., J. Phys. A., 46 (2013). [86] S. Muhuri and I. Pagonabarraga, Collective vesicle transport on biofilaments carried by competing molecular motors., EPL, 84 (2008). [87] M. J. I. Muller, S. Klumpp, and R. Lipowsky, Tug-of-war as a cooperative mechanism for bidirectional cargo transport by molecular motors., Proc. Natl. Acad. Sci. USA, 105 (2008), pp. 4609-4614. [88] J. M. Newby and P. C. Bressloff, Directed intermittent search for a hidden target on a dendritic tree., Phys. Rev. E, 80 (2009). [89] J. M. Newby and P. C. Bressloff, Quasi-steady state reduction of molecular-based models of directed intermittent search., Bull. Math. Biol., 72 (2010), pp. 1840-1866. [90] M. O'Toole, R. Latham, R. M. Baqri, and K. E. Miller, Modeling mitochondrial dynamics during in vivo elongation, J. Th. Biol., 255 (2008), pp. 369-377. [91] A. Parmeggiani, T. Franosch, and E. Frey, Totally asymmetric simple exclusion process with Langmuir kinetics, Phys. Rev. E, 70 (2004). [92] I. Pinkoviezky and N. S. Gov, Modeling interacting molecular motors with an internal degree of freedom, New J. Phys., 15 (2013). [93] I. Pinkoviezky and N. S. Gov, Transport dynamics of molecular motors that switch between an active and inactive state., Phys. Rev. E, 88 (2013). [94] A. E. Powell and A. Lenhard, Control of organ size in plants, Curr. Biol., 22 (2012), pp. R360-R367. [95] C. W. Pratt and K. Cornely, Essential Biochemistry, Wiley, 2010. [96] J. Prost, C. Barbetta, and J. F. Joanny, Dynamical control of the shape and size of stereocilia and microvilli., Biophys. J., 93 (2007), pp. 1124-1133. [97] S. M. Rafelski and W. F. Marshall, Building the cell: Design principles of cellular architecture, Mol. Cell Biol., 9 (2008), pp. 593-603. [98] M. C. Reed, S. Venakides, and J. J. Blum, Approximate traveling waves in linear reaction-hyperbolic equations., SIAM J. Appl. Math., 50 (1990), pp. 167-180. [99] L. Reese, A. Melbinger, and E. Frey, Crowding of molecular motors determines microtubule depolymerization., Biophys. J., 101 (2011), pp. 2190-2200. [100] T. Reichenbach, T. Franosch, and E. Frey, Exclusion processes with internal states., Phys. Rev. Lett., 97 (2006). [101] T. Reichenbach, E. Frey, and T. Franosch, Traffic jams induced by rare switching events in two-lane transport, New. J. Phys., 9 (2007). [102] I. Rishal, N. Kam, R. B. Perry, V. Shinder, E. M. C. Fisher, S. Giampietro, and M. Fainzilber, A motor-driven mechanism for cell-length sensing, Cell Reports, 1 (2012), pp. 608-616. 138 [103] A. K. Rzadzinska, M. E. Schneider, C. Davies, G. P. Riordan, and B. Kachar, An actin molecular treadmill and myosins maintain stereocilia functional architecture and self-renewal., J. Cell Biol., 64 (2004). [104] B. E. A. Saleh and M. C. Teich, Multiplied-Poisson noise in pulse, particle, and photon detection., Proc. IEEE, 70 (1982), pp. 229-245. [105] S. Salsa, Partial differential equations in action: from modelling to theory, Springer, 2015. [106] A. V. Samsonovich and G. A. Ascoli, Morphological homeostasis in cortical dendrites, Proc. Natl. Acad. Sci. USA, 103 (2006), pp. 1569-1574. [107] T. Sasamoto and H. Spohn, One-dimensional Kardar-Parisi-Zhang equation: an exact solution and its universality., Phys. Rev. Lett., 104.23 (2010). [108] J. M. Scholey, Intraflagellar transport., Ann. Rev. Cell Dev. Biol., 19 (2003), pp. 423- 443. [109] C. Settembre, A. Fraldi, D. L. Medina, and A. Ballabio, Signals from the lysosome: a control centre for cellular clearance and energy metabolism., Nature Reviews Molecular Cell Biology., 14 (2013), pp. 283-296. [110] S. Shinomoto and T. Yasuhiro, Modeling spiking behavior of neurons with timedependent Poisson processes., Phys. Rev. E, 64 (2001). [111] D. A. Smith and R. M. Simmons, Models of motor-assisted transport of intracellular particles., Biophys. J., 80 (2001), pp. 45-68. [112] D. H. Smith, Stretch growth of integrated axon tracts: extremes and exploitations, Prog. Neurobiol., 89 (2009), pp. 231-239. [113] R. E. Stephens, Quantal tektin synthesis and ciliary length in sea-urchin embryos., J. Cell Scie., 92 (1989), pp. 403-413. [114] C. Tischer, P. R. ten Wolde, and M. Dogterom, Providing positional information with active transport on dynamic microtubules., Biophys. J., 99 (2010), pp. 726-735. [115] M. Toriyama, Y. Sakumara, T. Shimada, S. Ishii, and N. Inagaki, A diffusion-based neurite length-sensing mechanism inolved in neuronal symmetry breaking., Molecular Systems Biology, 6 (2010). [116] F. Tostevin, Precision of sensing cell length via concentration gradients, Biophys J., 100 (2011), pp. 294-303. [117] F. Tostevin, W. de Ronde, and P. R. ten Wolde, Reliability of frequency and amplitude decoding in gene regulation, Phys. Rev. Lett., 108 (2012), p. 108104. [118] F. Tostevin, P. R. ten Wolde, and M. Howard, Fundamental limits to position determination by concentration gradients, PLoS Comput. Biol., 3 (2007), p. e78. [119] D. Tsygankov, A. W. R. Serohijos, N. V. Dokholyan, and T. C. Elston, A physical model reveals the mechanochemistry responsible for dynein's processive motion, Biophys J., 101 (2011), pp. 144-150. 139 [120] R. D. Vale, The molecular motor toolbox for intracellular transport., Cell., 112 (2003), pp. 467-480. [121] R. B. Vallee, J. C. Williams, D. Varma, and L. E. Barnhart, Dynein: An ancient motor protein involved in multiple modes of transport., Journal of Neurobiology, 58 (2004), pp. 189-200. [122] V. Varga, J. Helenius, K. Tanaka, A. A. Hyman, T. U. Tanaka, and J. Howard, Yeast kinesin-8 depolymerizes microtubules in a length-dependent manner., Nature Cell Biol., 8 (2006), pp. 957-962. [123] K. J. D. Vas, A. J. Grierson, S. Ackerly, and C. C. Miller, Role of axonal transport in neurodegenerative diseases., Ann. Rev. Neurosci., 31 (2008), pp. 151-173. [124] R. H. Wade, On and around microtubules: an overview., Mio. Biotechnol., 43 (2009), pp. 177-191. [125] M. Y. Wong, C. Zhou, D. Shakiryanova, T. E. Lloyd, D. L. Deitcher, and E. S. Levitan, Neuropeptide delivery to synapses by long-range vesicle circulation and sporadic capture., Cell, 148 (2012), pp. 1029-1038. [126] K. N. Wren, J. M. Craft, D. Tritschler, A. Schauer, D. K. Patel, E. F. Smith, M. E. Porter, P. Kner, and K. F. Lechtreck, A differential cargo-loading model of ciliary length regulation by IFT., Curr. Biol., 23 (2013), pp. 2463-2471. [127] S. Yogev, R. Cooper, R. Fetter, M. Horowitz, and K. Shen, Microtubule organization determines axonal transport dynamics, Neuron, 92 (2016), pp. 449-460. [128] S. Yukawa, M. Kikuchi, and S. Tadaki, Dynamical phase transition in one dimensional traffic flow model with blockage, Journal of the Physical Society of Japan, 63 (1994), pp. 3609-3618. SUBJECT INDEX active transport, 2, 16 Hill function, 97 adiabatic approximation, 10, 34, 57, 74, 96 Hopf bifurcation, 98 advection-diffusion equation, 17, 18, 23, 95, IFT, 71, 74 99 polar, 28 spherical, 31 Bessel function, 29 binomial distribution, 77 Cayley tree, 21 Chapman-Kolmogorov equation, 88, 108 kinesin, 4, 93 Lambert W function, 51 mean-field approximation, 49, 55, 95, 107 Michaelis-Menten, 97 microtubule, 2, 4 neuron, 2, 7 characteristic functional, 81 chemical master equation, 12, 76, 87, 118 continuum limit, 49, 55, 82 piecewise deterministic Markov process, 8, 95 Poisson process, 74 delay differential equation, 96 autocorrelation function, 76 dynein, 4, 93 Doubly stochastic, 80 Doubly stochastic solution, 80 exclusion process, 47, 52 power spectrum, 111 expectation, 48, 55, 81 synaptic democracy, 16, 19, 25, 29, 32, 35, 37, feed forward serial network, 115 Feynman-Kac formula, 88 63 system-size expansion, 77, 119 Fokker-Planck equation, 11, 34, 37, 78, 119 TASEP, 46, 59, 60 Gillespie algorithm, 63, 90 thinning algorithm, 85 Green's function, 101, 102 tug of war, 9 advection-diffusion equation, 125