A Computationally Eﬃcient Method for Solving Min System Reaction-Diﬀusion Equations within Growing and Dividing Domains that Approximate a Rod-Shaped Bacterium by William Christopher Carlquist A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in THE FACULTY OF GRADUATE STUDIES (Mathematics) The University of British Columbia (Vancouver) April 2012 c William Christopher Carlquist, 2012 � Abstract Few biological systems, showing complex pattern formation that spans multiple spatial and temporal scales, have been reduced in understanding to several components. The Min system in Escherichia coli, consists of three proteins, MinD, MinE, and MinC. Through ordered, cyclic, membrane binding and unbinding, facilitated by ATP hydrolysis, the Min system regulates the site of cell division in vivo. The Min system is tightly coupled with cell growth and division. Various mathematical models have been proposed to describe specific biological phenomena, arising from the Min system, but no model has been tested in a biologically relevant space, which grows and divides. In this thesis, I develop a computationally eﬃcient method for solving reaction-diﬀusion equations, in a growing and dividing geometry, which approximates growth and division in a rod-shaped bacterium. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Numerical Solutions to Reaction-Diﬀusion Equations in RodCell Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Previously Published Numerical Approaches for Solving Deterministic Min System Reaction-Diﬀusion Models in vivo . . . . . . . . . . 2.2 Outline of the Base 2 Multi-Time-Step Forward Euler Method and Description within this Chapter . . . . . . . . . . . . . . . . . . . . 2.3 Diﬀusion in Rod-Cell Geometry . . . . . . . . . . . . . . . . . . . . 2.4 Finite Diﬀerences in Rod-Cell Geometry . . . . . . . . . . . . . . . 2.4.1 Deriving Spatial Finite Diﬀerence Formulas . . . . . . . . . 2.4.2 Matching and Boundary Conditions . . . . . . . . . . . . . . 2.4.3 The base 2 Multi-Time-Step Forward Euler Method . . . . . 2.4.4 Computational Eﬃciency of the Base 2 Multi-Time-Step Forward Euler Method . . . . . . . . . . . . . . . . . . . . . . . 2.4.5 Conservation in the Base 2 Multi-Time-Step Forward Euler Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii 1 5 . 5 . 6 . 6 . 8 . 10 . 23 . 28 . 32 . 37 2.4.6 2.5 Spatial Convergence of the Multi-Time-Step Forward Euler Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.7 Temporal Convergence of the Multi-Time-Step Forward Euler Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.8 Convergence of the Multi-Time-Step Forward Euler Method to the Infinite domain solution . . . . . . . . . . . . . . . . . 2.4.9 Diﬀusion on the Rod-Cell Boundary . . . . . . . . . . . . . . Reactions in Rod-Cell Geometry . . . . . . . . . . . . . . . . . . . . . 39 . 42 . 43 . 45 . 48 3 Numerical Solutions to Reaction-Diﬀusion Equations on Changing Domains in Rod-Cell Geometry . . . . . . . . . . . . . . . . . 3.1 Previously Published Results and Numerical Methods for Solving Reaction-Diﬀusion Systems in Changing Domains . . . . . . . . . . . 3.2 Numerically Solving Reaction-Diﬀusion Equations in Growing RodCell Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Rod-Cell Growth by Adding Cylindrical Sections . . . . . . . 3.2.2 Rod-Cell Growth by Stretching the Rod-Cell Body . . . . . . 3.2.3 A Comparison of Rod-Cell Growing methods . . . . . . . . . . 3.2.4 Hybrid Growth, by Adding Cylindrical Sections and Stretching 3.3 Dividing Rod-Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 The Base 2 Multi-Time-Step Forward Euler Method for Solving Diﬀusion Equations in Dividing Rod-Cells . . . . . . . . . 3.3.2 The Base 2 Multi-Time-Step Forward Euler Method for Solving Diﬀusion Equations on Dividing Rod-Cell Boundaries . . . 3.4 Corrections for Changes in Total Amounts of Conserved Quantities during Rod-Cell Growth and Division . . . . . . . . . . . . . . . . . . 4 Examples of Reaction-Diﬀusion Solutions in Rod-Cell Geometry 4.1 The Huang et al. Model . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 The Huang et al. Model in a Wild Type Cell that Grows and Divides 4.3 The Huang et al. Model in a Continuously Growing Filamentous Mutant Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Future Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv 49 49 50 51 51 52 52 53 54 58 61 62 62 64 67 69 70 List of Tables 2.1 2.2 2.3 2.4 Relative errors for various time-stepping methods based on time step Base 2 multi-time-step forward Euler method sparsity. . . . . . . . . Computational time to calculate Matrices required for various methods. Computational time to calculate one numerical second using various methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v 34 35 36 36 List of Figures 2.1 2.2 2.3 Spatial convergence of the base 2 multi-time-step forward Euler method in rod-cell geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 Temporal convergence of the base 2 multi-time-step forward Euler method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Nonlinear fit of numerical solutions to the infinite domain solution . . 45 3.1 Grid point locations in r relative to the leading edge of the invaginated boundary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.1 4.2 Huang et al. solutions in a wild type cell that grows and divides . . . 65 The Huang et al. model in a continuously growing filamentous mutant cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 vi Acknowledgements Many thanks to my supervisor, Eric Cytrynbaum, for his clear guidance and encouragement, Dan Coombs for his thorough reading and thoughtful suggestions, and my family and girlfriend for much love and support. William Christopher Carlquist The University of British Columbia (Vancouver) April 2012 vii for Namy viii Chapter 1 Introduction Interconnection and complexity impedes biological understanding. Few spontaneous biological patterns, spanning multiple temporal and spatial scales, have been understandably reduced to several components. The Min system in Escherichia coli is one. The Min system, consisting of three known proteins, MinD, MinE, and MinC, regulates cell division by disrupting formation of the FtsZ ring, a contractile ring that divides a mother cell into two daughter cells, (de Boer et al., 1989). According to current understanding, MinD binds to ATP within the cytosol, forming MinD:ATP complexes, which dimerize (Szeto et al., 2003). Cytosolic MinD:ATP dimers then bind to the cell membrane, preferentially in the presence of membranebound MinD:ATP complexes, likely in polymeric or oligomeric form (Lackner et al., 2003), (Zonglin et al., 2002). MinC proteins, the system component responsible for FtsZ ring disruption, bind to membrane-bound MinD:ATP complexes (Lackner et al., 2003). MinE proteins also bind to membrane-bound MinD:ATP complexes, at a binding site that overlaps the MinC binding site (Luyan Ma and King, 2004), and may dislodge previously bound MinC proteins (Lackner et al., 2003). Once bound, MinE stimulates MinD ATPase activity. MinD dephosphorylates ATP into ADP, breaks membrane and MinE bonds, returning to the cytosol as a MinD:ADP complex (Zonglin et al., 2002). The cytosolic MinD:ADP complex then exchanges ADP for ATP, and the binding cycle continues. The Min protein binding cycle, in conjunction with the bacterial cell shape, results in Min protein stochastic poleto-pole switching, in cells shorter than roughly 2.5µm, and regular pole-to-pole oscillation, with a period of about 80s, in cells longer than roughly 3µm (FischerFriedrich et al., 2010), with maximal mean concentrations over time, at cell poles, and minimum mean concentrations over time, at midcell. At midcell, the minimum 1 mean MinC concentration over time results in a minimum disruption of FtsZ ring formation over time. Thus the FtsZ ring forms at midcell, dividing the mother cell into two roughly symmetric, viable daughter cells. In the absence of MinC, MinD and MinE oscillate in vivo. Also, the Min proteins decay slowly (Raskin and de Boer, 1999), so abrupt changes in protein concentration are not forcing oscillatory behavior. Additionally, cells lacking FstZ grow longer than wild type cells. At a certain length, roughly 7-8µm, Min oscillations bifurcate from a single pole-to-pole oscillation to two pole-to-midcell oscillations (Raskin and de Boer, 1999), indicating that a specific MinD binding partner, at cell poles, is likely not determining spatial organization in the Min system. Thus, the attempt to understand the oscillatory behavior of the Min system has been reduced to two proteins, MinD and MinE, an energy source, ATP, and the cell shape. In the Loose et al. in vitro experiments, fluorescently labeled MinD and MinE, with ATP, atop a supported lipid bilayer, produced MinD and MinE traveling and spiral waves (Loose et al., 2008), (Loose et al., 2011). Thus, complex spatial and temporal patterning emerges from only two proteins, MinD and MinE, an energy source, ATP, and a binding surface, the supported lipid bilayer. In a similar experiment, performed by Ivanov and Mizuuchi, a well mixed solution, containing ATP and fluorescently labeled MinD and MinE, was flowed atop a supported lipid bilayer. Dynamic patterns, such as traveling waves, spiral waves, propagating snake-like structures, and amoeba-like structures were observed (Ivanov and Mizuuchi, 2010). Complex spatial and temporal patterns, arising below the well mixed, constantly flowing solution, demonstrate that spatial and temporal variation in solution components are not driving Min protein pattern formation. Various mathematical models have been proposed to describe the Min system. Only one model, the Loose model, has attempted to describe the in vitro observations of the Loose et al. experiments (Loose et al., 2008). The model demonstrates traveling and spiral wave solutions, but traveling wave solutions diﬀer significantly in shape from experimental observations; shapes of spiral wave solutions have not been compared to experimental findings. Additionally, the model requires depletion of proteins within the buﬀer atop the lipid bilayer, a condition incompatible with the Ivanov and Mizuuchi findings. In a one dimensional in vitro context, the Loose model shows regular pole-to-pole oscillations for 2µm long cells, but other results have not been shown. Most mathematical models describing the Min system have focused on specific in vivo behaviors. The Cytrynbaum and Marshall model, a multistranded polymer 2 model, captures both regular pole-to-pole Min protein oscillations and the behaviors of various, experimentally tested, Min mutant phenotypes (Cytrynbaum and Marshall, 2007). Meinhardt and de Boer (Meinhardt and de Boer, 2001) modeled the Min system in a growing one dimensional cell. Their model demonstrates regular pole-to-pole oscillations that change into higher order oscillatory modes, during cell growth. However, their model requires fast protein production and degradation for oscillations, despite experiments showing persistent Min protein oscillations in the absence of Min protein production. The Huang et al model, a reaction diffusion model, shown in chapter 4, demonstrates regular pole-to-pole Min protein oscillations in 4µm long cells, and two regular pole-to-midcell Min protein oscillations in 10µm long cells (Huang et al., 2003), behavior consistent with experimentally demonstrated bifurcations in Min protein oscillation, occuring in long cells. A stochastic version of the Huang et al. model captures stochastic pole-to-pole Min protein switching in short cells and regular pole-to-pole Min protein oscillations in longer cells , behavior compatible with the Fisher-Friedrich et al. experimental findings (Fischer-Friedrich et al., 2010). Tostevin and Howard proposed a one dimensional stochastic MinD polymerization model, within a dividing domain, predicting that protein diﬀusion across the closing septum was unable to produce an equal distribution of Min proteins in daughter cells (Tostevin and Howard, 2006). In contrast, Ventura and Sourjik implemented a stochastic variation of the Huang et al. model, at various fixed septum closures, demonstrating a possible mechanism for equal Min protein distribution in daughter cells, a result matching their experimental findings (Ventura and Sourjik, 2011). Kerr et al. compared stochastic and deterministic versions of the Huang et al. model to study the Min system’s ability to find midcell (Kerr et al., 2006). Similarly, Fange and Elf compared stochastic and deterministic versions of the Huang et al model to address various mutant phenotypes (Fange and Elf, 2006). To study the Min system’s behavior during and after septation, Sengupta and Rutenburg deterministically solved the Huang et al. model in a dividing cell (Sengupta and Rutenburg, 2007). No model has addressed all in vitro and in vivo Min system behaviors. Numerically solving reaction-diﬀusion models in an in vitro setting is relatively straight forward. However, solving reaction-diﬀusion equations in a realistic in vivo setting, while maintaining computational eﬃciency, is more problematic. Only one model has addressed Min protein behavior in a three dimensional dividing cell, and no models have addressed Min protein behavior in a three dimensional growing cell. The Min system is tightly coupled with cell growth and division, thus to understand 3 model robustness, in a more realistic in vivo space, and to provide biologically pertinent predictions, model behaviors must be determined in three dimensional growing and dividing cells. This thesis outlines the computationally eﬃcient, base 2 multitime-step forward Euler method, for solving reaction-diﬀusion equations in growing and dividing rod-cell geometries, and provides examples and applications of numerical solutions found using the method. 4 Chapter 2 Numerical Solutions to Reaction-Diﬀusion Equations in Rod-Cell Geometry Eschericia Coli is a rod-shaped bacterium. To numerically solve reaction-diﬀusion equations, describing Min system protein dynamics, within an Eschericia Coli in vivo environment, I developed a computationally eﬃcient method, the base 2 multitime-step forward Euler method, within rod-cell geometry, a geometrical approximation to a rod-shaped bacterium. 2.1 Previously Published Numerical Approaches for Solving Deterministic Min System ReactionDiﬀusion Models in vivo Various geometries and numerical methods have been applied to solve deterministic reaction-diﬀusion models of the Min system, in vivo. Loose et al. (Loose et al., 2008) and Meinhardt and de Boer (Meinhardt and de Boer, 2001) modeled Eschericia Coli as a one dimensional line, and solved their systems of PDEs using the Forward Euler method. Huang et al. (Huang et al., 2003) and Sungupta and Rutenberg (Sengupta and Rutenburg, 2007) modeled Eschericia Coli as a cylinder, with equal grid point spacing in the r and in the z direction. They solved their PDE models using unmentioned time-stepping processes. Kerr et al. modeled Eschericia Coli 5 both as a box and as a cylinder, with equal grid point spacing throughout each geometry, and solved their model through simple explicit time-stepping (Kerr et al., 2006). Fange and Elf modeled Eschericia Coli by Constructive Solid Geometry, through the union of a hemispheres on each end of a cylinder (Fange and Elf, 2006), and numerically solved mean-field models through the method of lines, in conjunction with a two step backward diﬀerentiation formula. Like Fange and Elf, I model Eschericia Coli as the union of a hemisphere on each side of a cylinder, but provide a novel, computationally eﬃcient, time-stepping method for solving reaction-diﬀusion equations, the base 2 multi-time-step forward Euler method. 2.2 Outline of the Base 2 Multi-Time-Step Forward Euler Method and Description within this Chapter The base 2 multi-time-step forward Euler method for solving reaction-diﬀusion equations in rod-cell geometry implements multiple time steps, for computationally efficient solution calculations. To solve the diﬀusion parts of reaction-diﬀusion equations, detailed in section 2.4, the method utilizes finite diﬀerence approximation, described in section 2.4.1, with matching and boundary conditions, defined in section 2.4.2. Variable time-stepping, within the method, is described in section 2.4.3. Resulting computational eﬃciency discussion, and comparison of computational efficiency using other methods are shown in section 2.4.4. Solutions, found using the method, are shown to be conservative in section 2.4.5, spatially convergent in section 2.4.6, temporally convergent in section 2.4.7, and convergent to the known infinite domain solution, under similar conditions, in section 2.4.8. The method also extends to solving reaction diﬀusion equations on the rod-cell boundary, which is discussed in section 2.4.9. Finally, complete solutions to reaction-diﬀusion equations, in rod cell geometry, are calculated by adding changes from reactions to diﬀusion solutions, in section 2.5. 2.3 Diﬀusion in Rod-Cell Geometry Reaction-diﬀusion equations are of type: ∂f (x, t) = F (f (x, t)) + D∇2 f (x, t) , ∂t 6 (2.1) where f (x, t) represents concentrations of reaction species, at location x, at time t, F represents local reactions of reaction species, and D is a diagonal matrix of diﬀusion constants. To numerically solve systems of type Eq. (2.1), I split each equation within the system into reaction and diﬀusion parts, then developed a computationally eﬃcient method, the base 2 multi-time-step forward Euler method, to numerically solve the diﬀusion part of each equation within the system: ∂f (x, t) = D∇2 f (x, t) , ∂t (2.2) in rod-cell geometry. Rod-cell geometry consists of a hemisphere convexly adjoined to each side of a cylinder. To numerically solve the system of diﬀusion equations, Eq. (2.2), within rod-cell geometry, I separated the geometry into three natural parts—top and bottom caps, consisting of hemispheres, and the body, consisting of a cylinder. Separated into body and caps, the system of diﬀusion equations, Eq. (2.2), in rod-cell geometry, becomes: � 2 � ∂fi (r, θ, z, t) ∂ fi (r, θ, z, t) 1 ∂fi (r, θ, z, t) 1 ∂ 2 fi (r, θ, z, t) ∂ 2 fi (r, θ, z, t) = Di + + 2 + , ∂t ∂r2 r ∂r r ∂θ2 ∂z 2 (2.3) within the rod-cell body, � 2 ∂gi (r, θ, φ, t) ∂ gi (r, θ, φ, t) 2 ∂gi (r, θ, φ, t) 1 ∂ 2 gi (r, θ, φ, t) =Di + + ∂t ∂r2 r ∂r ∂θ2 r2 sin2 φ � cos φ ∂gi (r, θ, φ, t) 1 ∂ 2 gi (r, θ, φ, t) + 2 + 2 , (2.4) r sin φ ∂φ r ∂φ2 within the bottom cap, and � 2 ∂hi (r, θ, φ, t) ∂ hi (r, θ, φ, t) 2 ∂hi (r, θ, φ, t) 1 ∂ 2 hi (r, θ, φ, t) =Di + + ∂t ∂r2 r ∂r ∂θ2 r2 sin2 φ � cos φ ∂hi (r, θ, φ, t) 1 ∂ 2 hi (r, θ, φ, t) + 2 + 2 , (2.5) r sin φ ∂φ r ∂φ2 within the top cap, for each reaction species, i, of body length L, and radius R. Points at the interface of the body and either bottom or top cap are represented within both the body and the cap, hence: π fi (r, θ, 0) = gi (r, θ , ), 2 7 (2.6) at the bottom cap-body interface, and π fi (r, θ, L) = hi (r, θ , ), 2 (2.7) at the top cap-body interface. At the bottom cap-body interface, the first and second partial derivatives of f and g in z, represented, respectively, in cylindrical and spherical coordinates, should match, hence: ∂fi (r, θ, 0) 1 ∂gi (r, θ , π2 ) = , ∂z r ∂φ ∂ 2 fi (r, θ, 0) 1 ∂gi (r, θ , π2 ) 1 ∂ 2 gi (r, θ , π2 ) = + 2 , ∂z r ∂r r ∂ 2 φ2 (2.8) (2.9) Similarly, for the top cap-body interface: ∂fi (r, θ, L) 1 ∂hi (r, θ , π2 ) = , ∂z r ∂φ ∂ 2 fi (r, θ, L) 1 ∂hi (r, θ , π2 ) 1 ∂ 2 hi (r, θ , π2 ) , = + 2 ∂z r ∂r r ∂ 2 φ2 (2.10) (2.11) Fixing interface matching conditions, Eqs. 2.6, 2.8, 2.9, 2.10, 2.7, and 2.11 forces invariance in rod-cell diﬀusion equations, Eqs. 2.3 and 2.4, at the bottom cap-body interface, and Eqs. 2.3 and 2.5, at the top cap-body interface, resulting in smooth transitions across both cap-body interfaces. Reaction species are conserved within the rod-cell system, hence do not cross the rod-cell boundary. Boundary conditions are thus closed at r = R: ∂fi (R, θ, z, t) = 0, ∂r ∂gi (R, θ, φ, t) = 0, ∂r ∂hi (R, θ, φ, t) = 0. ∂r 2.4 (2.12) (2.13) (2.14) Finite Diﬀerences in Rod-Cell Geometry The system of diﬀusion equations, Eq. (2.2), within rod-cell geometry, separated into the body, bottom, and top cap, Eqs (2.3), (2.4), and (2.5), is numerically solved 8 by finite diﬀerence approximation. Finite diﬀerence approximation is chosen rather than finite element approximation, as rod-cell geometry is a union of simple geometries, which is uniformly discretized fairly easily, and finite diﬀerence implementation over various time-steps, for computational eﬃciency, is more straightforward than finite element implementation over various time-steps. Each spatial variable, r, θ, z, and φ, is uniformly spaced to define a three dimensional grid, on which finite diﬀerences are calculated. Within the rod-cell bottom cap, body, and top cap, there are Nr diﬀering grid points in r, and Nθ diﬀering grid points in θ. Within the rod-cell body, there are Nz diﬀering grid points in z. Within each cap, there are Nφ diﬀering grid points. A global grid density determining value, Nglobal , relates the number of grid points in each spatial dimension, such that: Nr = 2RNglobal , (2.15) Nθ = 4Nglobal , (2.16) Nz = 2LNglobal + 1, (2.17) Nφ = Nglobal + 1, (2.18) for rod-cell of radius, R ∈ {.5, 1, 1.5, 2, . . . }, body length, L ∈ {.5, 1, 1.5, 2, . . . }, and Nglobal ≥ 5. Grid points in r, of r-index j, are ordered from 1, at the r grid point nearest to r = 0, to Nr , at r = R, the rod-cell radius. In rod-cell geometry, singularities exist at r = 0 in the system of diﬀusion equations, Eqs. (2.3), (2.4), and (2.5). To avoid singular finite diﬀerence values and maintain equal grid point spacing across r = 0, the first grid point in r, of r-index 1, is located at r = ∆r , where ∆r is the 2 grid point spacing in r. Hence, ∆r + (Nr − 1)∆r = R. 2 � R ∆r = . Nr − 12 (2.19) (2.20) Grid points in θ, of θ-index k, are ordered from 1, at θ = 0, to Nθ , at θ = 2π − ∆θ, where ∆θ is the grid point spacing in θ, hence: ∆θ = 9 2π . Nθ (2.21) Grid points in z, of z-index l, are ordered from 1, at z = 0, to Nz , at z = L, hence: L . (2.22) Nz − 1 Within the bottom cap, grid points in φ, of φ-index m, are ordered from 1, at the grid point nearest to φ = 0, to Nφ , at φ = π2 . Within the top cap, grid points in φ, of φ-index m, are ordered from 1, at the grid point nearest to φ = π, to Nφ , at φ = π2 . In rod-cell geometry, singularities exist within the system of diﬀusion equations, Eqs. (2.3), (2.4), and (2.5), at φ = 0 and φ = π. To avoid singular finite diﬀerence values and maintain equal grid point spacing across φ = 0, in the bottom cap, and φ = π, in the top cap, the first grid point in φ, of φ-index 1, is located at φ = ∆φ , in the bottom cap, and φ = π − ∆φ , in the top cap, where ∆φ is the grid 2 2 point spacing in φ. Hence, ∆z = ∆φ π + (Nφ − 1)∆φ = . 2 2 � ∆φ = 2.4.1 π 2 Nφ − 1 2 . (2.23) (2.24) Deriving Spatial Finite Diﬀerence Formulas Changing variables in diﬀusion equations from cartesian to cylindrical or spherical geometries results in nonautonomous coeﬃcients, of the form 1r and r12 , in rod-cell cos φ 1 1 body diﬀusion equation, (2.3), and of the form 2r , r2 sin 2 φ , r 2 sin φ , and r 2 in rodcell cap diﬀusion equations, (2.4) and (2.5). These nonautonomous coeﬃcients can increase error in finite diﬀerence approximation unless handled carefully. Depending on the value of each nonautonomous coeﬃcient, finite diﬀerence approximations to spatial partial derivatives are calculated by scaling and combining Taylor expansions of the specific function, f, g, or h, in ψ, where ψ is any spatial variable, r, θ, z or φ, at the necessary number of evenly spaced grid points to attain a maximal finite diﬀerence approximation error of O(∆2 ), where ∆ = max{∆r, ∆θ, ∆z, ∆φ}, and O(∆r) = O(∆θ) = O(∆z) = O(∆φ). The second order partial derivative of f in z, in rod cell diﬀusion equation, (2.3), has a constant coeﬃcient, thus a symmetric three-point finite diﬀerence is 10 required for local finite diﬀerence approximation error to be of O(∆2 ): ∂ 2 f (r, θ, z, t) f (r, θ, z − ∆z, t) 2f (r, θ, z, t) f (r, θ, z + ∆z, t) ≈ − + ∂z 2 ∆z 2 ∆z 2 ∆z 2 2 4 ∆z ∂ f (r, θ, z, t) − . 12 ∂z 4 (2.25) Nonautonomous coeﬃcients for first order partial derivates of f , g, and h in r, in rod cell diﬀusion equations, (2.3), (2.4), and (2.5), are of O(r−1 ). At small r, nonautonomous coeﬃcients of O(r−1 ) increase error in finite diﬀerence approximation by one order of magnitude, O(∆), thus symmetric five point finite diﬀerences are required to approximate first order partial derivatives of f , g, and h, in r, with local finite diﬀerence approximation errors maximally of O(∆2 ). Finite diﬀerences that approximate second order partial derivatives of f , g, and h, in r, utilize the same grid points as finite diﬀerences approximating first order partial derivatives in r, hence: ∂f (r, θ, z, t) f (r − 2∆r, θ, z, t) 8f (r − ∆r, θ, z, t) 8f (r + ∆r, θ, z, t) ≈ − + ∂r 12∆r 12∆r 12∆r f (r + 2∆r, θ, z, t) ∆r4 ∂ 5 f (r, θ, z, t) − + , (2.26) 12∆r 30 ∂r5 ∂ 2 f (r, θ, z, t) −f (r − 2∆r, θ, z, t) 16f (r − ∆r, θ, z, t) 30f (r, θ, z, t) ≈ + − ∂r2 12∆r2 12∆r2 12∆r2 16f (r + ∆r, θ, z, t) f (r + 2∆r, θ, z, t) ∆r4 ∂ 5 f (r, θ, z, t) + − + , 12∆r2 12∆r2 90 ∂r6 (2.27) ∂g(r, θ, φ, t) g(r − 2∆r, θ, φ, t) 8g(r − ∆r, θ, φ, t) 8g(r + ∆r, θ, φ, t) ≈ − + ∂r 12∆r 12∆r 12∆r g(r + 2∆r, θ, φ, t) ∆r4 ∂ 5 g(r, θ, φ, t) − + , (2.28) 12∆r 30 ∂r5 ∂ 2 g(r, θ, φ, t) −g(r − 2∆r, θ, φ, t) 16g(r − ∆r, θ, φ, t) 30g(r, θ, φ, t) ≈ + − ∂r2 12∆r2 12∆r2 12∆r2 16g(r + ∆r, θ, φ, t) g(r + 2∆r, θ, φ, t) ∆r4 ∂ 5 g(r, θ, φ, t) + − + , 12∆r2 12∆r2 90 ∂r6 (2.29) 11 ∂h(r, θ, φ, t) h(r − 2∆r, θ, φ, t) 8h(r − ∆r, θ, φ, t) 8h(r + ∆r, θ, φ, t) ≈ − + ∂r 12∆r 12∆r 12∆r h(r + 2∆r, θ, φ, t) ∆r4 ∂ 5 h(r, θ, φ, t) − + , (2.30) 12∆r 30 ∂r5 ∂ 2 h(r, θ, φ, t) −h(r − 2∆r, θ, φ, t) 16h(r − ∆r, θ, φ, t) 30h(r, θ, φ, t) ≈ + − ∂r2 12∆r2 12∆r2 12∆r2 16h(r + ∆r, θ, φ, t) h(r + 2∆r, θ, φ, t) ∆r4 ∂ 5 h(r, θ, φ, t) + − + . 12∆r2 12∆r2 90 ∂r6 (2.31) Within the rod-cell body, grid points in r are located directly across the r origin, diﬀering by θ = π, thus may act as grid points in the negative r direction, for finite diﬀerence formulas in r, as done by Lai for a disk in polar coordinates (Lai, 2001). However, in both caps, there are no equivalent grid points located directly across the r origin, thus grid points in both caps of r index 1 and 2 must implement 5 point, non-symmetric finite diﬀerence formulas that require no grid points in the negative r direction. For r index 1: ∂g(r, θ, φ, t) −25g(r, θ, φ, t) 48g(r + ∆r, θ, φ, t) 36g(r + 2∆r, θ, φ, t) ≈ + − ∂r 12∆r 12∆r 12∆r 16g(r + 3∆r, θ, φ, t) 3g(r + 4∆r, θ, φ, t) ∆r4 ∂ 5 g(r, θ, φ, t) + − + , 12∆r 12∆r 5 ∂r5 (2.32) ∂ 2 g(r, θ, φ, t) 35g(r, θ, φ, t) 104g(r + ∆r, θ, φ, t) 114g(r + 2∆r, θ, φ, t) ≈ − + ∂r2 12∆r2 12∆r2 12∆r2 56g(r + 3∆r, θ, φ, t) 11g(r + 4∆r, θ, φ, t) 5∆r4 ∂ 5 g(r, θ, φ, t) − + + , 12∆r2 12∆r2 6 ∂r5 (2.33) ∂h(r, θ, φ, t) −25h(r, θ, φ, t) 48h(r + ∆r, θ, φ, t) 36h(r + 2∆r, θ, φ, t) ≈ + − ∂r 12∆r 12∆r 12∆r 16h(r + 3∆r, θ, φ, t) 3h(r + 4∆r, θ, φ, t) ∆r4 ∂ 5 h(r, θ, φ, t) + − + , 12∆r 12∆r 5 ∂r5 (2.34) ∂ 2 h(r, θ, φ, t) 35h(r, θ, φ, t) 104h(r + ∆r, θ, φ, t) 114h(r + 2∆r, θ, φ, t) ≈ − + ∂r2 12∆r2 12∆r2 12∆r2 56h(r + 3∆r, θ, φ, t) 11h(r + 4∆r, θ, φ, t) 5∆r4 ∂ 5 h(r, θ, φ, t) − + − . 12∆r2 12∆r2 6 ∂r5 (2.35) 12 For r index 2: ∂g(r, θ, φ, t) −3g(r − ∆r, θ, φ, t) 10g(r, θ, φ, t) 18g(r + ∆r, θ, φ, t) ≈ − + ∂r 12∆r 12∆r 12∆r 6g(r + 2∆r, θ, φ, t) g(r + 3∆r, θ, φ, t) ∆r4 ∂ 5 g(r, θ, φ, t) − + + , 12∆r 12∆r 20 ∂r5 (2.36) ∂ 2 g(r, θ, φ, t) 11g(r − ∆r, θ, φ, t) 20g(r, θ, φ, t) 6g(r + ∆r, θ, φ, t) ≈ − + ∂r2 12∆r2 12∆r2 12∆r2 4g(r + 2∆r, θ, φ, t) g(r + 3∆r, θ, φ, t) ∆r3 ∂ 5 g(r, θ, φ, t) + − − , 12∆r2 12∆r2 12 ∂r5 (2.37) ∂h(r, θ, φ, t) −3h(r − ∆r, θ, φ, t) 10h(r, θ, φ, t) 18h(r + ∆r, θ, φ, t) ≈ − + ∂r 12∆r 12∆r 12∆r 6h(r + 2∆r, θ, φ, t) h(r + 3∆r, θ, φ, t) ∆r4 ∂ 5 h(r, θ, φ, t) − + + , 12∆r 12∆r 20 ∂r5 (2.38) ∂ 2 h(r, θ, φ, t) 11h(r − ∆r, θ, φ, t) 20h(r, θ, φ, t) 6h(r + ∆r, θ, φ, t) ≈ − + ∂r2 12∆r2 12∆r2 12∆r2 4h(r + 2∆r, θ, φ, t) h(r + 3∆r, θ, φ, t) ∆r3 ∂ 5 h(r, θ, φ, t) + − − . 12∆r2 12∆r2 12 ∂r5 (2.39) The nonautonomous coeﬃcient for the second order partial derivative of f in θ, in rod cell diﬀusion equation, (2.3), is of O(r−2 ), which, for small r, increase local error from finite diﬀerence approximation by two orders of magnitude, O(∆2 ), thus a symmetric five-point finite diﬀerence is required to approximate the second order partial derivative of f in θ with local finite diﬀerence approximation error of O(∆2 ): ∂ 2 f (r, θ, z, t) −f (r, θ − 2∆θ, z, t) 16f (r, θ − ∆θ, z, t) 30f (r, θ, z, t) ≈ + − ∂θ2 12∆θ2 12∆θ2 12∆θ2 16f (r, θ + ∆θ, z, t) f (r, θ + 2∆θ, z, t) ∆θ4 ∂ 5 f (r, θ, z, t) + − + . 12∆θ2 12∆θ2 90 ∂θ6 (2.40) Nonautonomous coeﬃcients for second order partial derivative of g and h in θ, in rod cell diﬀusion equations (2.4) and (2.5), r−2 sin−2 (φ), increases local 13 error from finite diﬀerence approximation, at small r and small φ, by four orders of magnitude, O(∆4 ), thus symmetric seven-point finite diﬀerences are required to approximate second order partial derivatives of g and h in θ with local finite diﬀerence approximation errors maximally of O(∆2 ): ∂ 2 g(r, θ, φ, t) 2g(r, θ − 3∆θ, φ, t) 27g(r, θ − 2∆θ, φ, t) 270g(r, θ − ∆θ, φ, t) ≈ − + ∂θ2 180∆θ2 180∆θ2 180∆θ2 490g(r, θ, φ, t) 270g(r, θ + ∆θ, φ, t) 27g(r, θ + 2∆θ, φ, t) − + − 180∆θ2 180∆θ2 180∆θ2 6 8 2g(r, θ + 3∆θ, φ, t) ∆θ ∂ g(r, θ, φ, t) + − . (2.41) 180∆θ2 560 ∂θ8 ∂ 2 h(r, θ, φ, t) 2h(r, θ − 3∆θ, φ, t) 27h(r, θ − 2∆θ, φ, t) 270h(r, θ − ∆θ, φ, t) ≈ − + ∂θ2 180∆θ2 180∆θ2 180∆θ2 490h(r, θ, φ, t) 270h(r, θ + ∆θ, φ, t) 27h(r, θ + 2∆θ, φ, t) + − − 180∆θ2 180∆θ2 180∆θ2 6 8 2h(r, θ + 3∆θ, φ, t) ∆θ ∂ h(r, θ, φ, t) + − . (2.42) 180∆θ2 560 ∂θ8 Partial derivatives of g and h in φ, in rod cell diﬀusion equations (2.4) and (2.5), have nonautonomous coeﬃcients r2cos(φ) and r−2 . At small r and small φ, sin(φ) the nonautonomous coeﬃcient r2cos(φ) increases local errors from finite diﬀerence sin(φ) approximation by three orders of magnitude, thus symmetric seven-point finite differences of g and h in φ are required for local finite diﬀerence approximation errors to be maximally of O(∆2 ). Within each cap, grid points of φ-index less than Nφ − 2, have at least three additional φ grid points in the direction of the cap-body interface, thus symmetric seven-point finite diﬀerences that approximate first and second 14 order partial derivatives of g and h in φ are implemented: ∂g(r, θ, φ, t) −g(r, θ, φ − 3∆φ, t) 9g(r, θ, φ − 2∆φ, t) 45g(r, θ, φ − ∆φ, t) ≈ + − ∂φ 60∆φ2 60∆φ 60∆φ 45g(r, θ, φ + ∆φ, t) 9g(r, θ, φ + 2∆φ, t) g(r, θ, φ + 3∆φ, t) + − + 60∆φ 60∆φ 60∆φ ∆φ6 ∂ 7 g(r, θ, φ, t) − , (2.43) 140 ∂φ7 ∂ 2 g(r, θ, φ, t) 2g(r, θ, φ − 3∆φ, t) 27g(r, θ, φ − 2∆φ, t) 270g(r, θ, φ − ∆φ, t) ≈ − + ∂φ2 180∆φ2 180∆φ2 180∆φ2 490g(r, θ, φ, t) 270g(r, θ, φ + ∆φ, t) 27g(r, θ, φ + 2∆φ, t) − + − 180∆φ2 180∆φ2 180∆φ2 6 8 2g(r, θ, φ + 3∆φ, t) ∆φ ∂ g(r, θ, φ, t) + − , (2.44) 180∆φ2 560 ∂φ8 ∂h(r, θ, φ, t) −h(r, θ, φ − 3∆φ, t) 9h(r, θ, φ − 2∆φ, t) 45h(r, θ, φ − ∆φ, t) ≈ + − ∂φ 60∆φ2 60∆φ 60∆φ 45h(r, θ, φ + ∆φ, t) 9h(r, θ, φ + 2∆φ, t) h(r, θ, φ + 3∆φ, t) + − + 60∆φ 60∆φ 60∆φ 6 7 ∆φ ∂ h(r, θ, φ, t) − , (2.45) 140 ∂φ7 ∂ 2 h(r, θ, φ, t) 2h(r, θ, φ − 3∆φ, t) 27h(r, θ, φ − 2∆φ, t) 270h(r, θ, φ − ∆φ, t) ≈ − + ∂φ2 180∆φ2 180∆φ2 180∆φ2 490h(r, θ, φ, t) 270h(r, θ, φ + ∆φ, t) 27h(r, θ, φ + 2∆φ, t) − + − 180∆φ2 180∆φ2 180∆φ2 6 8 2h(r, θ, φ + 3∆φ, t) ∆φ ∂ h(r, θ, φ, t) + − , (2.46) 180∆φ2 560 ∂φ8 Within each cap, at grid points of φ-index, Nφ − 2, the nonautonomous coef−2 ficient, rcos(φ) ), thus symmetric five-point finite diﬀerences that approx2 sin(φ , is of O(r imate first and second order partial derivatives of g and h in φ, in rod cell diﬀusion equations (2.4) and (2.5), are required for local finite diﬀerence approximation errors 15 to be maximally of O(∆2 ): ∂g(r, θ, φ, t) g(r, θ, φ − 2∆φ, t) 8g(r, θ, φ − ∆φ, t) 8g(r, θ, φ + ∆φ, t) ≈ − + ∂φ 12∆φ 12∆φ 12∆φ 4 5 g(r, θ, φ + 2∆φ, t) ∆φ ∂ g(r, θ, φ, t) − + , (2.47) 12∆φ 30 ∂φ5 ∂ 2 g(r, θ, φ, t) −g(r, θ, φ − 2∆φ, t) 16g(r, θ, φ − ∆φ, t) 30g(r, θ, φ, t) ≈ + − ∂φ2 12∆φ2 12∆φ2 12∆φ2 16g(r, θ, φ + ∆φ, t) g(r, θ, φ + 2∆φ, t) ∆φ4 ∂ 6 g(r, θ, φ, t) + − + , 12∆φ2 12∆φ2 90 ∂φ6 (2.48) ∂h(r, θ, φ, t) h(r, θ, φ − 2∆φ, t) 8h(r, θ, φ − ∆φ, t) 8h(r, θ, φ + ∆φ, t) ≈ − + ∂φ 12∆φ 12∆φ 12∆φ 4 5 h(r, θ, φ + 2∆φ, t) ∆φ ∂ h(r, θ, φ, t) − + , (2.49) 12∆φ 30 ∂φ5 ∂ 2 h(r, θ, φ, t) −h(r, θ, φ − 2∆φ, t) 16h(r, θ, φ − ∆φ, t) 30h(r, θ, φ, t) ≈ + − ∂φ2 12∆φ2 12∆φ2 12∆φ2 16h(r, θ, φ + ∆φ, t) h(r, θ, φ + 2∆φ, t) ∆φ4 ∂ 6 h(r, θ, φ, t) + − + , 12∆φ2 12∆φ2 90 ∂φ6 (2.50) At each cap-body interface, two interface matching conditions, continuity of first and second order partial derivatives in z, allow the calculation of one finite diﬀerence function value across the interface, for finite diﬀerences in z and φ (derived in Section 2.4.2). Thus, finite diﬀerences in z and φ are calculated with at most one finite diﬀerence function value across the interface. Within each cap, at grid points of φ-index, Nφ − 1, the nonlinear partial derivative multiplier, r2cos(φ) , is of O(r−2 ). sin(φ) Implementing five point finite diﬀerences of g and h in φ, Eqs 2.47, 2.48, 2.49, and 2.50, forces local finite diﬀerence approximation errors to be maximally of O(∆2 ), and require only one finite diﬀerence function value across each interface. Within each cap, at grid points within the cap-body interface, of φ-index, Nφ , the nonautonomous coeﬃcients for the second order partial derivatives of g and h in φ, in rod cell diﬀusion equations (2.4) and (2.5), are of O(r−2 ). Finite diﬀerences may only implement one function value across each interface, calculated from interface matching conditions. Thus, non-symmetric six-point finite diﬀerences are required to approximate first and second order partial derivatives of g and h in 16 φ, with local finite diﬀerence approximation errors maximally of O(∆2 ): ∂g(r, θ, φ, t) 3g(r, θ, φ − 4∆φ, t) 20g(r, θ, φ − 3∆φ, t) 60g(r, θ, φ − 2∆φ, t) ≈ − + ∂φ 60∆φ2 60∆φ2 60∆φ 120g(r, θ, φ − ∆φ, t) 65g(r, θ, φ, t) 12g(r, θ, φ + ∆φ, t) − + + 60∆φ 60∆φ 60∆φ 5 6 ∆φ ∂ g(r, θ, φ, t) − , (2.51) 30 ∂φ6 ∂ 2 g(r, θ, φ, t) g(r, θ, φ − 4∆φ, t) 6g(r, θ, φ − 3∆φ, t) 14g(r, θ, φ − 2∆φ, t) ≈ − + ∂φ2 12∆φ2 12∆φ2 12∆φ2 4g(r, θ, φ − ∆φ, t) 15g(r, θ, φ, t) 10g(r, θ, φ + ∆φ, t) − − + 12∆φ2 12∆φ2 12∆φ2 4 6 13∆φ ∂ g(r, θ, φ, t) − , (2.52) 180 ∂φ6 ∂h(r, θ, φ, t) 3h(r, θ, φ − 4∆φ, t) 20h(r, θ, φ − 3∆φ, t) 60h(r, θ, φ − 2∆φ, t) ≈ − + ∂φ 60∆φ2 60∆φ2 60∆φ 120h(r, θ, φ − ∆φ, t) 65h(r, θ, φ, t) 12h(r, θ, φ + ∆φ, t) − + + 60∆φ 60∆φ 60∆φ 5 6 ∆φ ∂ h(r, θ, φ, t) − , (2.53) 30 ∂φ6 ∂ 2 h(r, θ, φ, t) h(r, θ, φ − 4∆φ, t) 6h(r, θ, φ − 3∆φ, t) 14h(r, θ, φ − 2∆φ, t) ≈ − + ∂φ2 12∆φ2 12∆φ2 12∆φ2 4h(r, θ, φ − ∆φ, t) 15h(r, θ, φ, t) 10h(r, θ, φ + ∆φ, t) − − + 12∆φ2 12∆φ2 12∆φ2 4 6 13∆φ ∂ h(r, θ, φ, t) − , (2.54) 180 ∂φ6 Finite Diﬀerence Definitions For reaction species i, at grid points of r-index j, θ-index k, z-index l, φ-index m, and temporal, t,-index n, the following are defined: ∆2r fi (rj , θk , zl , tn ) = −fi (rj−2 , θk , zl , tn ) 16fi (rj−1 , θk , zl , tn ) 30fi (rj , θk , zl , tn ) + − 12∆r2 12∆r2 12∆r2 16fi (rj+1 , θk , zl , tn ) fi (rj+2 , θk , zl , tn ) + − , (2.55) 12∆r2 12∆r2 17 for j ∈ {1, . . . , Nr }. ∆1r fi (rj , θk , zl , tn ) = for j ∈ {1, . . . , Nr }. ∆2θ fi (rj , θk , zl , tn ) = fi (rj−2 , θk , zl , tn ) 8fi (rj−1 , θk , zl , tn ) 8fi (rj+1 , θk , zl , tn ) − + 12∆r 12∆r 12∆r (2.56) fi (rj+2 , θk , zl , tn ) − , (2.57) 12∆r fi (rj , θk−2 , zl , tn ) 16fi (rj , θk−1 , zl , tn ) 30fi (rj , θk , zl , tn ) + − 12∆θ2 12∆θ2 12∆θ2 16fi (rj , θk+1 , zl , tn ) fi (rj , θk+2 , zl , tn ) + − , (2.58) 12∆θ2 12∆θ2 for k ∈ {1, . . . , Nθ }. ∆2z fi (rj , θk , zl , tn ) = for l ∈ {1, . . . , Nz }. ∆2r gi (rj , θk , φm , tn ) = for j ∈ {1}. ∆2r gi (rj , θk , φm , tn ) = for j ∈ {2}. ∆2r gi (rj , θk , φm , tn ) = for j ∈ {3, . . . , Nr }. ∆1r gi (rj , θk , φm , tn ) = fi (rj , θk , zl−1 , tn ) 2fi (rj , θk , zl , tn ) fi (rj , θk , zl+1 , tn ) − + , ∆z 2 ∆z 2 ∆z 2 (2.59) 35gi (rj , θk , φm , tn ) 104gi (rj+1 , θk , φm , tn ) 114gi (rj+2 , θk , φm , tn ) − + 12∆r2 12∆r2 12∆r2 56gi (rj+3 , θk , φm , tn ) 11gi (rj+4 , θk , φm , tn ) − + , (2.60) 12∆r2 12∆r2 11gi (rj−1 , θk , φm , tn ) 20gi (rj , θk , φm , tn ) 6gi (rj+1 , θk , φm , tn ) − + 12∆r2 12∆r2 12∆r2 4gi (rj+2 , θk , φm , tn ) gi (rj+3 , θk , φm , tn ) + − , (2.61) 12∆r2 12∆r2 −gi (rj−2 , θk , φm , tn ) 16gi (rj−1 , θk , φm , tn ) 30gi (rj , θk , φm , tn ) + − 12∆r2 12∆r2 12∆r2 16gi (rj+1 , θk , φm , tn ) gi (rj+2 , θk , φm , tn ) + − , (2.62) 12∆r2 12∆r2 −25gi (rj , θk , φm , tn ) 48gi (rj+1 , θk , φm , tn ) 36gi (rj+2 , θk , φm , tn ) + − 12∆r 12∆r 12∆r 16gi (rj+3 , θk , φm , tn ) 3gi (rj+4 , θk , φm , tn ) + − , (2.63) 12∆r 12∆r 18 for j ∈ {1}. ∆1r gi (rj , θk , φm , tn ) = for j ∈ {2}. ∆1r gi (rj , θk , φm , tn ) = for j ∈ {3, . . . , Nr }. ∆2θ gi (rj , θk , φm , tn ) = for k ∈ {1, . . . , Nθ }. ∆1φ gi (rj , θk , φm , tn ) = −3gi (rj−1 , θk , φm , tn ) 10gi (rj , θk , φm , tn ) 18gi (rj+1 , θk , φm , tn ) − + 12∆r 12∆r 12∆r 6gi (rj+2 , θk , φm , tn ) gi (rj+3 , θk , φm , tn ) − + , (2.64) 12∆r 12∆r gi (rj−2 , θk , φm , tn ) 8gi (rj−1 , θk , φm , tn ) 8gi (rj+1 , θk , φm , tn ) − + 12∆r 12∆r 12∆r (2.65) gi (rj+2 , θk , φm , tn ) − , (2.66) 12∆r 2gi (rj , θk−3 , φm , tn ) 27gi (rj , θk−2 , φm , tn ) 270gi (rj , θk−1 , φm , tn ) − + 180∆θ2 180∆θ2 180∆θ2 490gi (rj , θk , φm , tn ) 270gi (rj , θk+1 , φm , tn ) 27gi (rj , θk+2 , φm , tn ) − + − 180∆θ2 180∆θ2 180∆θ2 2gi (rj , θk+3 , φm , tn ) + , (2.67) 180∆θ2 −gi (rj , θk , φm−3 , tn ) 9gi (rj , θk , φm−2 , tn ) 45gi (rj , θk , φm−1 , tn ) + − 60∆φ 60∆φ 60∆φ 45gi (rj , θk , φm+1 , tn ) 9gi (rj , θk , φm+2 , tn ) gi (rj , θk , φm+3 , tn ) − − + , 60∆φ 60∆φ 60∆φ (2.68) for m ∈ {1, . . . , Nφ − 3}. ∆1φ gi (rj , θk , φm , tn ) = gi (rj , θk , φm−2 , tn ) 8gi (rj , θk , φm−1 , tn ) 8gi (rj , θk , φm+1 , tn ) − + 12∆φ 12∆φ 12∆φ gi (rj , θk , φm+2 , tn ) − , (2.69) 12∆φ for m ∈ {Nφ − 2, Nφ − 1}. ∆1φ gi (rj , θk , φm , tn ) = 3gi (rj , θk , φm−4 , tn ) 20gi (rj , θk , φm−3 , tn ) 60gi (rj , θk , φm−2 , tn ) − + 60∆φ 60∆φ 60∆φ 120gi (rj , θk , φm−1 , tn ) 65gi (rj , θk , φm , tn ) 12gi (rj , θk , φm+1 , tn ) − + + , 60∆φ 60∆φ 60∆φ (2.70) 19 for m ∈ {Nφ }. ∆2φ gi (rj , θk , φm , tn ) = 2gi (rj , θk , φm−3 , tn ) 27gi (rj , θk , φm−2 , tn ) 270gi (rj , θk , φm−1 , tn ) − + 180∆φ2 180∆φ2 180∆φ2 490gi (rj , θk , φm , tn ) 270gi (rj , θk , φm+1 , tn ) 27gi (rj , θk , φm+2 , tn ) − + − 180∆φ2 180∆φ2 180∆φ2 2gi (rj , θk , φm+3 , tn ) + , (2.71) 180∆φ2 for m ∈ {1, . . . , Nφ − 3}. ∆2φ gi (rj , θk , φm , tn ) = −gi (rj , θk , φm−2 , tn ) 16gi (rj , θk , φm−1 , tn ) 30gi (rj , θk , φm , tn ) + − 12∆φ2 12∆φ2 12∆φ2 16gi (rj , θk , φm+1 , tn ) gi (rj , θk , φm+2 , tn ) + − , (2.72) 12∆φ2 12∆φ2 for m ∈ {Nφ − 2, Nφ − 1}. ∆2φ gi (rj , θk , φm , tn ) = gi (rj , θk , φm−4 , tn ) 6gi (rj , θk , φm−3 , tn ) 14gi (rj , θk , φm−2 , tn ) − + 12∆φ2 12∆φ2 12∆φ2 4gi (rj , θk , φm−1 , tn ) 15gi (rj , θk , φm , tn ) 10gi (rj , θk , φm+1 , tn ) − − + , 12∆φ2 12∆φ2 12∆φ2 (2.73) for m ∈ {Nφ }. ∆2r hi (rj , θk , φm , tn ) = 35hi (rj , θk , φm , tn ) 104hi (rj+1 , θk , φm , tn ) 114hi (rj+2 , θk , φm , tn ) − + 12∆r2 12∆r2 12∆r2 56hi (rj+3 , θk , φm , tn ) 11hi (rj+4 , θk , φm , tn ) − + , (2.74) 12∆r2 12∆r2 for j ∈ {1}. ∆2r hi (rj , θk , φm , tn ) = 11hi (rj−1 , θk , φm , tn ) 20hi (rj , θk , φm , tn ) 6hi (rj+1 , θk , φm , tn ) − + 12∆r2 12∆r2 12∆r2 4hi (rj+2 , θk , φm , tn ) hi (rj+3 , θk , φm , tn ) + − , (2.75) 12∆r2 12∆r2 for j ∈ {2}. ∆2r hi (rj , θk , φm , tn ) = −hi (rj−2 , θk , φm , tn ) 16hi (rj−1 , θk , φm , tn ) 30hi (rj , θk , φm , tn ) + − 12∆r2 12∆r2 12∆r2 16hi (rj+1 , θk , φm , tn ) hi (rj+2 , θk , φm , tn ) + − , (2.76) 12∆r2 12∆r2 20 for j ∈ {3, . . . , Nr }. ∆1r hi (rj , θk , φm , tn ) = −25hi (rj , θk , φm , tn ) 48hi (rj+1 , θk , φm , tn ) 36hi (rj+2 , θk , φm , tn ) + − 12∆r 12∆r 12∆r 16hi (rj+3 , θk , φm , tn ) 3hi (rj+4 , θk , φm , tn ) + − , (2.77) 12∆r 12∆r for j ∈ {1}. ∆1r hi (rj , θk , φm , tn ) = −3hi (rj−1 , θk , φm , tn ) 10hi (rj , θk , φm , tn ) 18hi (rj+1 , θk , φm , tn ) − + 12∆r 12∆r 12∆r 6hi (rj+2 , θk , φm , tn ) hi (rj+3 , θk , φm , tn ) − + , (2.78) 12∆r 12∆r for j ∈ {2}. ∆1r hi (rj , θk , φm , tn ) = hi (rj−2 , θk , φm , tn ) 8hi (rj−1 , θk , φm , tn ) 8hi (rj+1 , θk , φm , tn ) − + 12∆r 12∆r 12∆r (2.79) hi (rj+2 , θk , φm , tn ) − , (2.80) 12∆r for j ∈ {3, . . . , Nr }. ∆2θ hi (rj , θk , φm , tn ) = 2hi (rj , θk−3 , φm , tn ) 27hi (rj , θk−2 , φm , tn ) 270hi (rj , θk−1 , φm , tn ) − + 180∆θ2 180∆θ2 180∆θ2 490hi (rj , θk , φm , tn ) 270hi (rj , θk+1 , φm , tn ) 27hi (rj , θk+2 , φm , tn ) − + − 180∆θ2 180∆θ2 180∆θ2 2hi (rj , θk+3 , φm , tn ) + , (2.81) 180∆θ2 for k ∈ {1, . . . , Nθ }. ∆1φ hi (rj , θk , φm , tn ) = −hi (rj , θk , φm−3 , tn ) 9hi (rj , θk , φm−2 , tn ) 45hi (rj , θk , φm−1 , tn ) + − 60∆φ 60∆φ 60∆φ 45hi (rj , θk , φm+1 , tn ) 9hi (rj , θk , φm+2 , tn ) hi (rj , θk , φm+3 , tn ) − − + , 60∆φ 60∆φ 60∆φ (2.82) for m ∈ {1, . . . , Nφ − 3}. ∆1φ hi (rj , θk , φm , tn ) = hi (rj , θk , φm−2 , tn ) 8hi (rj , θk , φm−1 , tn ) 8hi (rj , θk , φm+1 , tn ) − + 12∆φ 12∆φ 12∆φ hi (rj , θk , φm+2 , tn ) − , (2.83) 12∆φ 21 for m ∈ {Nφ − 2, Nφ − 1}. ∆1φ hi (rj , θk , φm , tn ) = 3hi (rj , θk , φm−4 , tn ) 20hi (rj , θk , φm−3 , tn ) 60hi (rj , θk , φm−2 , tn ) − + 60∆φ 60∆φ 60∆φ 120hi (rj , θk , φm−1 , tn ) 65hi (rj , θk , φm , tn ) 12hi (rj , θk , φm+1 , tn ) − + + , 60∆φ 60∆φ 60∆φ (2.84) for m ∈ {Nφ }. ∆2φ hi (rj , θk , φm , tn ) = 2hi (rj , θk , φm−3 , tn ) 27hi (rj , θk , φm−2 , tn ) 270hi (rj , θk , φm−1 , tn ) − + 180∆φ2 180∆φ2 180∆φ2 490hi (rj , θk , φm , tn ) 270hi (rj , θk , φm+1 , tn ) 27hi (rj , θk , φm+2 , tn ) − + − 180∆φ2 180∆φ2 180∆φ2 2hi (rj , θk , φm+3 , tn ) + , (2.85) 180∆φ2 for m ∈ {1, . . . , Nφ − 3}. ∆2φ hi (rj , θk , φm , tn ) = −hi (rj , θk , φm−2 , tn ) 16hi (rj , θk , φm−1 , tn ) 30hi (rj , θk , φm , tn ) + − 12∆φ2 12∆φ2 12∆φ2 16hi (rj , θk , φm+1 , tn ) hi (rj , θk , φm+2 , tn ) + − , (2.86) 12∆φ2 12∆φ2 for m ∈ {Nφ − 2, Nφ − 1}. ∆2φ hi (rj , θk , φm , tn ) = hi (rj , θk , φm−4 , tn ) 6hi (rj , θk , φm−3 , tn ) 14hi (rj , θk , φm−2 , tn ) − + 12∆φ2 12∆φ2 12∆φ2 4hi (rj , θk , φm−1 , tn ) 15hi (rj , θk , φm , tn ) 10hi (rj , θk , φm+1 , tn ) − − + , 12∆φ2 12∆φ2 12∆φ2 (2.87) for m ∈ {Nφ }. Spatial Finite Diﬀerences Grid points at the interface of the rod-cell body and either cap are represented in both cylindrical and spherical coordinates, thus diﬀusion occurs within both cylindrical and spherical geometry, and is invariant from matching conditions. To 22 ensure that diﬀusivity is not double the intended diﬀusivity, the total diﬀusion, at the interface of the rod-cell body and either cap, is half cylindrical diﬀusion and half spherical diﬀusion. Spatial finite diﬀerence formulas for the system of diﬀusion equations, Eq. (2.2), in rod-cell geometry, divided into the rod-cell body, Eq. (2.3), bottom cap, Eq. (2.4), and top cap, Eq. (2.5), for reaction species i, at grid points of r-index j, θ-index k, z-index l, φ-index m, and temporal, t,-index n are thus: ∂fi (rj , θk , zl , tn ) ≈ ∂t �� � � δl1 + δlNz 1 Di 1 − · ∆2r fi (rj , θk , zl , tn ) + ∆1r fi (rj , θk , zl , tn ) 2 rj �� 1 2 2 + 2 ∆θ fi (rj , θk , zl , tn ) + ∆z fi (rj , θk , zl , tn ) rj (2.88) ∂gi (rj , θk , φm , tn ) ≈ ∂t � � �� δmNφ 2 Di 1 − · ∆2r gi (rj , θk , φm , tn ) + ∆1r gi (rj , θk , φm , tn ) 2 rj 1 cos φm 1 + 2 2 ∆2θ gi (rj , θk , φm , tn ) + 2 ∆1φ gi (rj , θk , φm , tn ) + 2 ∆2φ gi (rj , θk , φm , tn ) rj sin φm rj rj sin φm �� (2.89) ∂hi (rj , θk , φm , tn ) ≈ ∂t �� � � δmNφ 2 Di 1 − · ∆2r hi (rj , θk , φm , tn ) + ∆1r hi (rj , θk , φm , tn ) 2 rj 1 cos φm 1 + 2 2 ∆2θ hi (rj , θk , φm , tn ) + 2 ∆1φ hi (rj , θk , φm , tn ) + 2 ∆2φ hi (rj , θk , φm , tn ) rj sin φm rj rj sin φm (2.90) where δxy is the Kronecker delta, with spatial variable indices x and y. 2.4.2 Matching and Boundary Conditions Finite diﬀerence boundary conditions are chosen to match those described in section 2.3, for diﬀusion in rod-cell geometry. Grid points at cap-body interfaces are 23 �� represented in both the body and the cap, hence: fi (rj , θk , z1 ) = gi (rj , θj , φNφ ), (2.91) at the bottom cap-body interface, and fi (rj , θk , z1 ) = hi (rj , θj , φNφ ), (2.92) FInite diﬀerence approximations of bottom cap-body matching conditions, 2.8 and 2.9, matching first and second order partial derivative in z, implementing the appropriate number of grid points for approximations to be maximally of O(∆2 ): 0= fi (rj , θk , z0 , tn ) fi (rj , θk , z2 , tn ) 3gi (rj , θk , φNφ −4 , tn ) 20gi (rj , θk , φNφ −3 , tn ) − + − 2∆z 2∆z 60rj ∆φ 60rj ∆φ 60gi (rj , θk , φNφ −2 , tn ) 120gi (rj , θk , φNφ −1 , tn ) 65gi (rj , θk , φNφ , tn ) + − + 60rj ∆φ 60rj ∆φ 60rj ∆φ 12gi (rj , θk , φNφ +1 , tn ) + , (2.93) 60rj ∆φ and 0= −fi (rj , θk , z0 , tn ) 2fi (rj , θk , z1 , tn ) fi (rj , θk , z2 , tn ) gi (rj−2 , θk , φNφ , tn ) + − + ∆z 2 ∆z 2 ∆z 2 12rj ∆r 8gi (rj−1 , θk , φNφ , tn ) 8gi (rj+1 , θk , φNφ , tn ) gi (rj+2 , θk , φNφ , tn ) − + − 12rj ∆r 12rj ∆r 12rj ∆r gi (rj , θk , φNφ −4 , tn ) 6gi (rj , θk , φNφ −3 , tn ) 14gi (rj , θk , φNφ −2 , tn ) + − + 12rj2 ∆φ2 12rj2 ∆φ2 12rj2 ∆φ2 4gi (rj , θk , φNφ −1 , tn ) 15gi (rj , θk , φNφ , tn ) 10gi (rj , θk , φNφ +1 , tn ) − − + , (2.94) 12rj2 ∆φ2 12rj2 ∆φ2 12rj2 ∆φ2 are scaled and combined to define function values that cross the bottom cap-body interface, fi (rj , θk , z0 , tn ) and gi (rj , θk , φNφ +1 , tn ), in terms of known function values, 24 which are then implemented in finite diﬀerence formulas. � � 24rj ∆φ 50∆z fi (rj ,θk , z0 , tn ) = fi (rj , θk , z1 , tn ) + + 1 fi (rj , θk , z2 , tn ) 25∆z + 12rj ∆φ 25∆z + 12rj ∆φ 3∆z 2 32∆z 2 − gi (rj , θk , φNφ −4 , tn ) + gi (rj , θk , φNφ −3 , tn ) 2rj (25∆z + 12rj ∆φ) 3rj (25∆z + 12rj ∆φ) 36∆z 2 96∆z 2 − gi (rj , θk , φNφ −2 , tn ) + gi (rj , θk , φNφ −1 , tn ) rj (25∆z + 12rj ∆φ) rj (25∆z + 12rj ∆φ) 415∆z 2 ∆φ∆z 2 − gi (rj , θk , φNφ , tn ) + gi (rj−2 , θk , φNφ , tn ) 6rj (25∆z + 12rj ∆φ) ∆r(25∆z + 12rj ∆φ) 8∆φ∆z 2 8∆φ∆z 2 − gi (rj−1 , θk , φNφ , tn ) + gi (rj+1 , θk , φNφ , tn ) ∆r(25∆z + 12rj ∆φ) ∆r(25∆z + 12rj ∆φ) ∆φ∆z 2 − gi (rj+2 , θk , φNφ , tn ) (2.95) ∆r(25∆z + 12rj ∆φ) −60rj2 ∆φ2 gi (rj , θk , φNφ +1 , tn ) = fi (rj , θk , z1 , tn ) ∆z(25∆z + 12rj ∆φ) � � 60rj2 ∆φ2 15∆z 1 gi (rj , θk , φNφ −4 , tn ) + fi (rj , θk , z2 , tn ) + − ∆z(25∆z + 12rj ∆φ) 4(25∆z + 12rj ∆φ) 4 � � � � 80∆z 5 90∆z − + gi (rj , θk , φNφ −3 , tn ) + − 5 gi (rj , θk , φNφ −2 , tn ) 3(25∆z + 12rj ∆φ) 3 25∆z + 12rj ∆φ � � � � 240∆z 2075∆z 65 − + 10 gi (rj , θk , φNφ −1 , tn ) + − gi (rj , θk , φNφ , tn ) 25∆z + 12rj ∆φ 12(25∆z + 12rj ∆φ) 12 5rj ∆φ2 ∆z 20rj ∆φ2 ∆z − gi (rj−2 , θk , φNφ , tn ) + gi (rj−1 , θk , φNφ , tn ) 2∆r(25∆z + 12rj ∆φ) ∆r(25∆z + 12rj ∆φ) 20rj ∆φ2 ∆z 5rj ∆φ2 ∆z − gi (rj+1 , θk , φNφ , tn ) + gi (rj+2 , θk , φNφ , tn ) ∆r(25∆z + 12rj ∆φ) 2∆r(25∆z + 12rj ∆φ) (2.96) Similarly, by combining finite diﬀerence approximations to matching conditions at the top cap-body interface, function values that cross the interface are defined in 25 terms of known function values: 24rj ∆φ fi (rj ,θk , zNz +1 , tn ) = fi (rj , θk , zNz , tn ) 25∆z + 12rj ∆φ � � 50∆z 3∆z 2 + + 1 fi (rj , θk , zNz −1 , tn ) − hi (rj , θk , φNφ −4 , tn ) 25∆z + 12rj ∆φ 2rj (25∆z + 12rj ∆φ) 32∆z 2 36∆z 2 + hi (rj , θk , φNφ −3 , tn ) − hi (rj , θk , φNφ −2 , tn ) 3rj (25∆z + 12rj ∆φ) rj (25∆z + 12rj ∆φ) 96∆z 2 415∆z 2 + hi (rj , θk , φNφ −1 , tn ) − hi (rj , θk , φNφ , tn ) rj (25∆z + 12rj ∆φ) 6rj (25∆z + 12rj ∆φ) ∆φ∆z 2 8∆φ∆z 2 + hi (rj−2 , θk , φNφ , tn ) − hi (rj−1 , θk , φNφ , tn ) ∆r(25∆z + 12rj ∆φ) ∆r(25∆z + 12rj ∆φ) ∆φ∆z 2 8∆φ∆z 2 hi (rj+1 , θk , φNφ , tn ) − hi (rj+2 , θk , φNφ , tn ) + ∆r(25∆z + 12rj ∆φ) ∆r(25∆z + 12rj ∆φ) (2.97) −60rj2 ∆φ2 hi (rj , θk , φNφ +1 , tn ) = fi (rj , θk , zNz , tn ) ∆z(25∆z + 12rj ∆φ) � � 60rj2 ∆φ2 15∆z 1 + fi (rj , θk , zNz −1 , tn ) + − hi (rj , θk , φNφ −4 , tn ) ∆z(25∆z + 12rj ∆φ) 4(25∆z + 12rj ∆φ) 4 � � � � 80∆z 5 90∆z + hi (rj , θk , φNφ −3 , tn ) + − 5 hi (rj , θk , φNφ −2 , tn ) − 3(25∆z + 12rj ∆φ) 3 25∆z + 12rj ∆φ � � � � 240∆z 2075∆z 65 − + 10 hi (rj , θk , φNφ −1 , tn ) + − hi (rj , θk , φNφ , tn ) 25∆z + 12rj ∆φ 12(25∆z + 12rj ∆φ) 12 5rj ∆φ2 ∆z 20rj ∆φ2 ∆z − hi (rj−2 , θk , φNφ , tn ) + hi (rj−1 , θk , φNφ , tn ) 2∆r(25∆z + 12rj ∆φ) ∆r(25∆z + 12rj ∆φ) 20rj ∆φ2 ∆z 5rj ∆φ2 ∆z − hi (rj+1 , θk , φNφ , tn ) + hi (rj+2 , θk , φNφ , tn ) ∆r(25∆z + 12rj ∆φ) 2∆r(25∆z + 12rj ∆φ) (2.98) Reaction species do not cross the rod-cell boundary, hence boundaries are closed at r = R, where r-index = Nr . To impose close boundaries at the rod-cell boundary, ghost points are introduced where r > R, at indices Nr+1 and Nr+2 . Fixing: fi (rNr +1 , θk , zl ) = fi (rNr , θk , zl ), (2.99) fi (rNr +2 , θk , zl ) = fi (rNr , θk , zl ), (2.100) 26 for each reaction species, i, within the rod-cell body, fixing: gi (rNr +1 , θk , φm ) = gi (rNr , θk , φm ), (2.101) gi (rNr +2 , θk , φm ) = gi (rNr , θk , φm ), (2.102) for each reaction species, i, within the bottom cap, and fixing: hi (rNr +1 , θk , φm ) = hi (rNr , θk , φm ), (2.103) hi (rNr +2 , θk , φm ) = hi (rNr , θk , φm ) (2.104) for each reaction species, i, within the top cap, imposes closed boundaries across r = R. ∆θ is chosen such that there are an even number of intervals dividing θ. Even division of θ ensures that all grid points at θ-index k have a corresponding grid point of opposite spatial orientation, diﬀering by θ = π, at θ-index (k + N2θ ) mod Nθ . This allows demarkation: fi (r0 , θk , zl ) = fi (r1 , θ(k+ Nθ ) mod N , zl ), (2.105) fi (r−1 , θk , zl ) = fi (r2 , θ(k+ Nθ ) mod N , zl ) (2.106) θ 2 θ 2 for each reaction species, i, within the rod-cell body. Similar to demarking function values at grid points of r-index 0 and −1, function values at grid points of φ-index 0, −1, and −2 are demarked: gi (rj , θk , φ0 ) = gi (rj , θ(k+ Nθ ) mod N , φ1 ), (2.107) gi (rj , θk , φ−1 ) = gi (rj , θ(k+ Nθ ) mod N , φ2 ), (2.108) gi (rj , θk , φ−2 ) = gi (rj , θ(k+ Nθ ) mod N , φ3 ), (2.109) θ 2 θ 2 θ 2 for each reaction species, i, within the bottom cap, and hi (rj , θk , φ0 ) = hi (rj , θ(k+ Nθ ) mod N , φ1 ), (2.110) hi (rj , θk , φ−1 ) = hi (rj , θ(k+ Nθ ) mod N , φ2 ), (2.111) hi (rj , θk , φ−2 ) = hi (rj , θ(k+ Nθ ) mod N , φ3 ), (2.112) θ 2 2 2 27 θ θ for each reaction species, i, within the top cap. Grid points are periodic in θ, hence: fi (rj , θ0 , zl ) = fi (rj , θNθ , zl ), (2.113) fi (rj , θ−1 , zl ) = fi (rj , θNθ −1 , zl ), (2.114) fi (rj , θNθ +1 , zl ) = fi (rj , θ1 , zl ), (2.115) fi (rj , θNθ +2 , zl ) = fi (rj , θ2 , zl ), (2.116) for each reaction species, i, within the rod-cell body, gi (rj , θ0 , φm ) = gi (rj , θNθ , φm ), (2.117) gi (rj , θ−1 , φm ) = gi (rj , θNθ −1 , φm ), (2.118) gi (rj , θ−2 , φm ) = gi (rj , θNθ −2 , φm ), (2.119) gi (rj , θNθ +1 , φm ) = gi (rj , θ1 , φm ), (2.120) gi (rj , θNθ +2 , φm ) = gi (rj , θ2 , φm ), (2.121) gi (rj , θNθ +3 , φm ) = gi (rj , θ3 , φm ), (2.122) for each reaction species, i, within the bottom cap, and hi (rj , θ0 , φm ) = hi (rj , θNθ , φm ), (2.123) hi (rj , θ−1 , φm ) = hi (rj , θNθ −1 , φm ), (2.124) hi (rj , θ−2 , φm ) = hi (rj , θNθ −2 , φm ), (2.125) hi (rj , θNθ +1 , φm ) = hi (rj , θ1 , φm ), (2.126) hi (rj , θNθ +2 , φm ) = hi (rj , θ2 , φm ), (2.127) hi (rj , θNθ +3 , φm ) = hi (rj , θ3 , φm ), (2.128) for each reaction species, i, within the top cap. 2.4.3 The base 2 Multi-Time-Step Forward Euler Method Spatial finite diﬀerences for diﬀusion in rod-cell geometry contain coeﬃcients that vary in orders of magnitude, based on spatial dimension and grid point location, arising from nonautonomous coeﬃcients in rod cell diﬀusion equations, (2.3), (2.4), and (2.5). Large coeﬃcient variation means that diﬀusion in rod-cell geometry evolves on diﬀerent time scales. The explicit forward Euler method can solve diﬀusion equations in rod cell geometry, but only stably when using a very short time 28 step, as large coeﬃcient variation results in instability under longer time steps. Implicit solution methods, such as the backward Euler method and Crank-Nicolson method, can numerically solve diﬀusion equations in rod cell geometry using longer time steps than the forward Euler method, but are computationally expensive. Naturally, an explicit numerical method that evolves only within pertinent time scales, on individual spatial dimensions, at individual grid points, would allow larger timestepping than the forward Euler method, while maintaining more computational eﬃciency than implicit methods. This method, called the base 2 multi-time-step forward Euler method, chooses variable time steps, ∆τs , for each spatial dimension, ψ ∈ {r, θ, z, φ}, at each grid point, such that the maximum of each spatial dimension finite diﬀerence coeﬃcient, cψ,i,j,k,l , in the rod-cell body, and cψ,i,j,k,m , within bottom and top caps, multiplied by the variable time step, ∆τs , are of equal order in magnitude, for reaction species i, at grid points of r-index j, θ-index k, z-index l, and φ index m. Additionally, at cap-body interfaces, finite diﬀerences in z evolve under the same variable time step as finite diﬀerences in φ for uniform time stepping throughout matching conditions. Within the rod-cell body: cr,i,j,k,l = Di · max{1, rj−1 }, cθ,i,j,k,l = Di · cz,i,j,k,l = � rj−2 , Di · rj−2 for l = 1 or l = Nz ; Di otherwise (2.129) (2.130) (2.131) Within both caps: cr,i,j,k,m = Di · max{1, 2rj−1 }, (2.132) cφ,i,j,k,m (2.134) cθ,i,j,k,m = Di · rj−2 2 sin φm , � � cos φm −2 = Di · max 2 ,r . rj sin φm j (2.133) The base 2 multi-time-step forward Euler method employs Nt diﬀerent time steps: ∆τs = ∆t 2Nt −s , (2.135) for s ∈ {1, . . . , Nt }, which cover the base 2 multi-time-step forward Euler macro time step, ∆t. 29 Finite diﬀerences, at each each grid point, in each spatial variable, evolve in only one time-step, chosen such that: Nt −s ≥ αi · cψ,i,j,k,l > 2Nt −s−1 ; 1 if 2 χ∆τs (cψ,i,j,k,l ) = 1 if αi · cψ,i,j,k,l ≤ 2−1 ; 0 otherwise, (2.136) within the rod-cell body, where αi is a fine tuning parameter that can shift finite diﬀerence evolutions under longer or shorter time steps, chosen for each reaction species, i, to ensure stability and maximum computational eﬃciency. Time steps, in both caps, at each grid point, for each spatial variable, are chosen such that: Nt −s ≥ αi · cψ,i,j,k,m > 2Nt −s−1 ; 1 if 2 χ∆τs (cψ,i,j,k,m ) = 1 if αi · cψ,i,j,k,m ≤ 2−1 ; 0 otherwise. (2.137) Rod-cell geometry base 2 forward Euler multi-time-step finite diﬀerences, for each time step ∆τs , of reaction species, i, including all grid points in the rod-cell bottom cap, body, and top cap, are thus : (fi (rj , θk , zl , tn + ∆τs ) + gi (rj , θk , φm , tn + ∆τs ) + hi (rj , θk , φm , tn + ∆τs )) − (fi (rj , θk , zl , tn ) + gi (rj , θk , φm , tn ) + hi (rj , θk , φm , tn )) = �� � � δl1 + δlNz ∆τs · Di 1 − · χ∆τs (Di · max{1, rj−1 }) · ∆2r fi (rj , θk , zl , tn ) 2 1 + χ∆τs (Di · max{1, rj−1 }) · ∆1r fi (rj , θk , zl , tn ) rj 1 + χ∆τs (Di · rj−2 ) · 2 ∆2θ fi (rj , θk , zl , tn ) rj + (δl1 + δlNz ) · χ∆τs (Di · rj−2 ) · ∆2z fi (rj , θk , zl , tn ) �� 2 + (1 − δl1 − δlNz ) · χ∆τs (Di ) · ∆z fi (rj , θk , zl , tn ) 30 +∆τs · Di �� δmNφ 1− 2 � � · χ∆τs (Di · max{1, 2rj−1 }) · ∆2r gi (rj , θk , φm , tn ) 2 1 ∆ gi (rj , θk , φm , tn ) rj r 1 + χ∆τs (Di · rj−2 sin−2 φm ) · 2 2 ∆2θ gi (rj , θk , φm , tn ) rj sin φm cos φm cos φm + χ∆τs (Di · max{ 2 , rj−2 }) · 2 ∆1 gi (rj , θk , φm , tn ) rj sin φm rj sin φm φ �� 1 2 cos φm −2 , r }) · 2 ∆φ gi (rj , θk , φm , tn ) + χ∆τs (Di · max{ 2 rj sin φm j rj �� � � δmNφ +∆τs · Di 1 − · χ∆τs (Di · max{1, 2rj−1 }) · ∆2r hi (rj , θk , φm , tn ) 2 2 + χ∆τs (Di · max{1, 2rj−1 }) · ∆1r hi (rj , θk , φm , tn ) rj 1 + χ∆τs (Di · rj−2 sin−2 φm ) · 2 2 ∆2θ hi (rj , θk , φm , tn ) rj sin φm cos φm cos φm + χ∆τs (Di · max{ 2 , rj−2 }) · 2 ∆1 hi (rj , θk , φm , tn ) rj sin φm rj sin φm φ �� cos φm 1 2 −2 + χ∆τs (Di · max{ 2 , r }) · 2 ∆φ hi (rj , θk , φm , tn ) rj sin φm j rj + χ∆τs (Di · max{1, 2rj−1 }) · (2.138) There are Nr Nθ (2Nφ + Nz − 2) = Ngrid total grid points in rod-cell geometry, ordered from 1 to Ngrid . Fi,s is the Ngrid × Ngrid sparse matrix, such that: Fi,s fi,t = ∆s fi,t , (2.139) where fi,t , is a column vector with Ngrid entries that correspond to the concentration of reaction species i, at time t, at each grid point, ordered from 1 to Ngrid , and ∆s fi,t is a column vector with Ngrid entries that correspond to the rod-cell geometry base 2 multi-time-step forward Euler finite diﬀerence value: (fi (rj , θk , zl , tn + ∆τs ) + gi (rj , θk , φm , tn + ∆τs ) + hi (rj , θk , φm , tn + ∆τs )) − (fi (rj , θk , zl , tn ) + gi (rj , θk , φm , tn ) + hi (rj , θk , φm , tn )) , for the time step, ∆τs , of reaction species, i, at time, t, at each grid point, ordered from 1 to Ngrid . 31 Complete finite diﬀerence evolutions are calculated by initially incorporating only the finite diﬀerence evolutions that occur during the shortest time step, ∆t1 , keeping constant all grid point values unaﬀected by finite diﬀerence evolutions during ∆t1 : fi,t+∆t1 = (I + Fi,1 )fi,t , (2.140) where I is a Ngrid × Ngrid sparse identity matrix. Time steps are chosen under a base 2 interval subdivision, thus finite diﬀerence evolutions under the shortest time step must occur a second time to cover the time interval, ∆t2 , the time interval under which the second set of finite diﬀerence evolutions occur. Thus, fi,t+∆t2 = ((I + Fi,1 )(I + Fi,1 ) + Fi,2 ) fi,t . Similarly, for all finite diﬀerence evolutions covering ∆t3 : �� �� � � fi,t+∆t3 = (I + Fi,1 )2 + Fi,2 (I + Fi,1 )2 + Fi,2 + Fi,3 fi,t . (2.141) (2.142) Continuing the iterative process, covering all finite diﬀerence evolutions under the base 2 multi-time-step forward Euler macro time step, ∆t: � �2 � � �2 � 2 � � 2 fi,t+∆t = . . . (I + Fi,1 )2 + Fi,2 + Fi,3 . . . + Fi,Nt−1 + Fi,Nt fi,t . (2.143) Ai , the complete base 2 multi-time-step forward Euler matrix, covering all finite diﬀerence evolutions over the macro time step ∆t, satisfies: fi,t+∆t = Ai fi,t , (2.144) and is calculated: � �2 � � �2 � 2 � � 2 Ai = . . . (I + Fi,1 )2 + Fi,2 + Fi,3 . . . + Fi,Nt−1 + Fi,Nt , (2.145) where I is a Ngrid × Ngrid sparse identity matrix. 2.4.4 Computational Eﬃciency of the Base 2 Multi-TimeStep Forward Euler Method Incorporating finite diﬀerence values only within pertinent time scales, while utilizing sparse matrix multiplication, prolongs sparsity in Ai . Calculating Ai under a 32 base n time step interval sub-division, where n > 1, would result in Ai calculations of the form: Ai = ��� n n . . . (((I + Fi,1 )n + Fi,2 ) + Fi,3 ) . . . �n + Fi,Nt−1 �n � + Fi,Nt . (2.146) Choosing time step intervals under a base 2 sub-division results in only one matrix multiplication before adding finite diﬀerences within the next time interval, as opposed to n − 1 matrix multiplications under a base n time step interval sub-division, where n > 1. Matrix multiplication of sparse Fi,s matrices reduces sparsity in the product matrix, and increases Ai computational time. A base 2 time step interval subdivision minimizes the number of matrix multiplications required to calculate Ai , maximizing sparsity in Ai , while minimizing Ai computational time, under matrix multiplication. The base 2 multi-time-step forward Euler macro time step, ∆t, and fine tuning parameter, αi , are chosen to optimize computational eﬃciency, based on the specific problem. For any ∆t, a αi can be chosen so that numerical solutions are stable and convergent. However, values of αi > 1 require longer Ai computational time and reduce Ai sparsity. For problems requiring few Ai calculations, over very extended time intervals, large values of ∆t coupled with large values of αi may be computationally beneficial, but for problems requiring many Ai calculations, such as rod-cell division or growing domain problems, computational time is minimized by choosing αi = 1 and finding the largest value of ∆t that returns stable, convergent solutions. Cell division and growing domain problems are the main foci of this thesis, thus quicker and sparse Ai computation, calculated under αi = 1, will be emphasized. To compare computational time requirements across various time-stepping methods, first, appropriate time steps for each method must be determined. For each time-stepping method, forward Euler, Backward Euler, Crank-Nicolson, and base 2 multi-time-step forward Euler, I calculated numerical solutions under time steps (macro time steps for the base 2 multi-time-step forward Euler method): 10−8 s, 10−7 s,. . . , 10−2 s. Then, I calculated the l 2 relative error for each solution compared to a solution found using the same method, under a time step one order of magnitude shorter, calculated over the same time interval, and starting from the same initial condition–a concentration value of 1000µm−3 at one grid point in the bottom cap, of r index 1, θ index 1, and φ index 1, with zero concentration values at all other grid points. All solutions were calculated in a rod-cell of length 1.5µm, under Di = 1µm2 s−1 , Nglobal = 5, and αi = 1. Results are shown in table 2.1. The maximum appropriate time step for each time-stepping method was chosen to 33 be the maximum time step with relative error less than .1. Thus for each method, backward Euler, Crank-Nicolson, forward Euler, and base 2 multi-time-step forward Euler, maximum appropriate time steps are, respectively, 10−6 s, 10−6 s, 10−7 s, and 10−3 s. Table 2.1: Relative errors for various time-stepping methods based on time step. Backward Euler is abbreviated BE. Crank-Nicolson is abbreviated CN. Forward Euler is abbreviated FE. Base 2 multi-time-step forward Euler is abbreviated MTS. A dash indicates unstable solutions. All solutions were calculated in a rod-cell of length 1.5µm, under Di = 1µm2 s−1 , Nglobal = 5, and αi = 1. Method time step BE CN FE MTS 10−8 s 10−7 s 10−6 s 10−5 s 10−4 s 10−3 s 10−2 s 2.1 · 10−5 1.9 · 10−3 .09 .34 .68 1.6 5.0 3.5 · 10−8 3.4 · 10−5 .02 - 2.1 · 10−5 2.2 · 10−3 .24 - 3.5 · 10−9 3.4 · 10−7 2.1 · 10−5 1.4 · 10−4 1.6 · 10−3 .02 .27 The base 2 multi-time-step forward Euler method maintains computational eﬃciency by evolving over a relatively long macro time step, while maintaing sparsity. For example, when Di = 1µm2 s−1 , Nglobal = 5, and αi = 1, the base 2 multi-time-step forward Euler method’s maximum appropriate macro time step is 10−3 s, while the forward Euler method’s maximum appropriate time step is 10−7 s. Additionally, the complete base 2 multi-time-step forward Euler finite diﬀerence matrix, Ai , contains only about an order of magnitude more non-zero matrix elements than the analogous forward Euler finite diﬀerence matrix, shown in Table 2.2 for various rod-cell lengths. Hence, the base 2 multi-time-step forward Euler method is approximately three orders of magnitude more computationally eﬃcient, when calculating diﬀusion solutions in rod-cell geometry, than the Forward Euler Method, shown in Table 2.4 for various rod-cell lengths, while requiring computational times of roughly the same order to calculate finite diﬀerence matrices, shown in Table 2.3 for various cell lengths. All computational times in this section were calculated in 34 MATLAB R2010a, on a MacBook Pro 7.1. Table 2.2: Base 2 multi-time-step forward Euler method sparsity. Sparsity of forward Euler and base 2 multi-time-step forward Euler finite diﬀerence matrices are shown. Forward Euler is abbreviated as FE. Base 2 multi-time-step forward Euler is abbreviated as MTS. Non-zero matrix elements is abbreviated as NZME. Calculations are from rod-cells with Di = 1µm2 s−1 , under a grid with Nglobal = 5, αi = 1, and base 2 multi-time-step forward Euler macro time step ∆t = 10−3 s. 2µm Total Matrix Elements 4.4 · 106 FE NZME. 2.9 · 104 MTS NZME. 9.2 · 105 FE fraction NZME. 6.5 · 10−3 MTS fraction NZME. 21.0 · 10−2 4µm Rod-cell length 6µm 8µm 16.8 · 106 4.9 · 104 11.0 · 105 2.9 · 10−3 6.5 · 10−2 37.2 · 106 7.0 · 104 12.7 · 105 1.9 · 10−3 3.4 · 10−2 65.6 · 106 9.0 · 104 14.5 · 105 1.4 · 10−3 2.2 · 10−2 10µm 102.0 · 106 11.0 · 104 16.2 · 105 1.1 · 10−3 1.6 · 10−2 The base 2 multi-time-step forward Euler method for solving diﬀusion equations in rod-cell geometry is computationally more eﬃcient than backward Euler or Crank-Nicolson implicit methods, as it implements sparse matrix multiplication rather than more computationally expensive implicit matrix solving methods, while covering a larger time interval per solution calculation. For example, when Di = 1µm2 s−1 and Nglobal = 5, the backward Euler and Crank-Nicolson methods’ maximum appropriate time step is 10−6 s. Comparatively, under the same grid, the base 2 multi-time-step forward Euler method’s maximum appropriate macro time step is 10−3 s, calculated under αi = 1, three orders of magnitude longer than the maximum appropriate time step for the backward Euler or Crank-Nicolson methods. Backward Euler and Crank-Nicolson implicit solutions can be found in two diﬀerent ways. Either solutions are found by solving a set of linear equations at each time step, or by generating a completely non-sparse matrix to calculate time step evolutions by matrix multiplication, computed by matrix inversion. Calculating general backward Euler and Crank-Nicolson matrices to solve diﬀusion equations by matrix multiplication is computationally expensive, but reduces solution computational time compared to solving a set of linear equations at each time step, for certain grids. The computational time to calculate matrices required for each 35 method, when Di = 1µm2 s−1 , Nglobal = 5, and αi = 1, are shown in Table 2.3, for various rod-cell lengths. The computational time for each method to calculate solutions up to 1s of numerical time, when Di = 1µm2 s−1 , Nglobal = 5, and αi = 1, are shown in Table 2.4, for rod-cells of various lengths. Table 2.3: Computational time to calculate Matrices required for various methods. BE Implicit and CN Implicit are backward Euler and Crank-Nicolson methods carrying out a linear solve at each time step. BEM and CNM are backward Euler and Crank-Nicolson methods where each solution is calculated by matrix multiplication. FE is the forward Euler method. MTS is the base 2 multi-time-step forward Euler method. Methods implement maximum appropriate time steps. Di = 1µm2 s−1 , Nglobal = 5, and αi = 1. Method 2µm Rod-cell length 4µm 6µm 8µm BE Implicit BEM CN Implicit CNM FE MTS .6s 9.9s .6s 9.7s .6s 5.7s .9s 1.3s 1.6s 2.1s 33.5s 139.9s 596.5s 864.7s .9s 1.3s 1.6s 2.1s 35.0s 143.4s 614.7s 892.2s .9s 1.3s 1.6s 2.1s 7.1s 8.2s 9.8s 11.1s 10µm Table 2.4: Computational time to calculate one numerical second using various methods. Definitions and parameter values are as in Table 2.3 Method 2µm BE Implicit 1.6 · 105 s BEM 1.3 · 104 s . CN Implicit 1.6 · 105 s CNM 1.3 · 104 s FE 7.0 · 102 s MTS 2.3s Rod-cell length 4µm 6µm 8µm 2.5 · 105 s 5.0 · 104 s 2.5 · 105 s 5.8 · 104 s 1.2 · 103 s 3.3s 8.4 · 105 s 1.1 · 105 s 9.1 · 105 s 1.1 · 105 s 1.7 · 103 s 3.9s 1.6 · 106 s 3.3 · 105 s 2.1 · 106 s 3.1 · 105 s 3.1 · 103 s 4.4s 10µm 3.0 · 106 s 6.0 · 105 s 2.9 · 106 s 5.9 · 105 s 3.8 · 103 s 5.0s Relative computational time amongst methods, as shown in Tables 2.3 and 2.4, is similar for grids with larger values of Nglobal and larger diﬀusion coeﬃcients, 36 Di . Of the methods shown, the base 2 multi-time-step forward Euler method is clearly the most computationally eﬃcient method for calculating diﬀusion in rodcells of various length. The biological time scale for cell division in Eschericia Coli is about twenty minutes. The base 2 multi-time-step forward Euler method makes computation on biologically relevant time scales tractable. 2.4.5 Conservation in the Base 2 Multi-Time-Step Forward Euler Method Diﬀusion in rod-cell geometry is a closed system. The total amount of each reaction species must be conserved during each time step of numerical diﬀusion. Let v be a scalar multiple of the vector of grid point volumes, and let fi,t be an arbitrary vector of concentrations, of reaction species i, at time t. Then, the total mass of reaction species i, at time t, is vT fi,t . The concentration of reaction species i, at time t + ∆t is Ai fi,t , with total mass vT Ai fi,t . The base 2 multi-time-step forward Euler Method is conserving, for each reaction species i, if and only if: vT Ai fi,t = vT fi,t , , (2.147) for any fi,t . vT Ai fi,t = vT fi,t � vT Ai fi,t − vT fi,t = 0 (2.148) � vT (Ai fi,t − fi,t ) = 0 (2.149) � vT (Ai − I)fi,t = 0 (2.150) � T fi,t (Ai − I)T v = 0 37 (2.151) � T fi,t (ATi − I T )v = 0 (2.152) � T fi,t (ATi − I)v = 0. (2.153) If 1 is an eigenvalue of ATi with corresponding eigenvector, v, then (ATi − I)v = 0. (2.154) � T fi,t (ATi − I)v = 0, (2.155) � vT Ai fi,t = vT fi,t for any fi,t . Thus, if ATi has eigenvalue, 1, and strictly positive or strictly negative (indicating a possible vector of concentrations) corresponding eigenvector, v, then the total amount of reaction species i is conserved during each time step of numerical diﬀusion, along v. The eigenvalue, 1, of ATi and corresponding strictly positive or strictly negative eigenvalue, v, are verified using the MATLAB function “eigs”, which numerically calculates eigenvalues of sparse matrices. If ATi has unique eigenvalue, 1, which is determined using the MATLAB function “eigs”, with corresponding strictly positive or strictly negative eigenvalue, v, then only scalings of v satisfy the conservation criterium, Eq. (2.147). An ordered column vector with Ngrid entries, of which each entry is the rod-cell numeric volume at a unique grid point, must satisfy the conservation criterium, Eq. (2.147), hence is a scaling of v. Therefore, v can be scaled into vvolume , an ordered column vector with Ngrid entries, of which each entry is the rod-cell numeric volume at a unique grid point. Hence, if vp is the pth entry of v: vvolume = rod-cell volume v, �Ngrid v p p=1 38 (2.156) and satisfies: T T vvolume Ai fi,t = vvolume fi,t , (2.157) vT fi,t = total amount of reaction species, i, and (2.158) Ngrid � vvolume,p = rod-cell volume, (2.159) p=1 for any fi,t , at all time, t, where vvolume,p is the pth entry of vvolume . 2.4.6 Spatial Convergence of the Multi-Time-Step Forward Euler Method Numerical solutions to the system of diﬀusion equations, Equation 2.2, in rod-cell geometry, using the base 2 multi-time-step forward Euler method, should converge spatially to an exact solution as the numeric-volume of each grid point decreases to zero. To test for spatial convergence, numerical solution were compared in six diﬀerent grids, (Nglobal = {5, 6, 7, 8, 9, 10}). There are four grid points common to all six grids: x1 (r-index = Nr , θ-index = 1, z-index = 1), x2 (r-index = Nr , θ-index = 1 + N2θ , z-index = 1), x3 (r-index = Nr , θ-index = 1, z-index = Nz ), and x4 (r-index = Nr , θ-index = 1 + N2θ , z-index = Nz ). A non-negative test function in cartesian coordinates, varying in all spatial directions, with diﬀering values at points x1 , x2 , x3 , and x4 : ftest = 100 + x + 50 cos(2π · x) cos(2π · y) cos(2π · z), (2.160) was used as the initial condition during the spatial convergence study. In the test function, the z origin is located where the bottom cap meets the body, at grid points of z-index 1. Spatial convergence was studied in rod-cells of body length, L = .5µm, with radius, R = .5µm, diﬀusion coeﬃcient, Di = 1µm2 s−1 , and under the time step fine-tuning parameter, αi = 1. The ordered vector containing the numericNglobal volume of each grid point, under Nglobal grid points, is denoted as vvolume . For each Nglobal Nglobal ∈ {5, 6, 7, 8, 9, 10}, vvolume · ftest = 91.629, where ftest is an ordered vector, of which each entry has value ftest , at a specific grid point location in cartesian coordinates. Hence, all initial conditions under grids with Nglobal ∈ {5, 6, 7, 8, 9, 10} have roughly equal total mass. Thus, there is little error from diﬀering total masses, when measuring spatial convergence. 10 n Numerical solution concentration diﬀerences, |fi,t (xj ) − fi,t (xj )|, 39 n ∈ {5, 6, 7, 8, 9}, were calculated at each grid point, x1 , x2 , x3 , and x4 , after one macro time step, ∆t = 10−4 s. Numerical solutions diﬀer at each grid point, x1 , x2 , x3 and x4 , but numerical concentration diﬀerences are identical, up to numerical roundoﬀ, at grid points x1 and x4 , and at grid points x2 and x3 . Symmetry within grids results in equal numeric-volume of all grid points, x1 , x2 , x3 , and x4 , within each grid selection. Figure 2.1 shows numerical concentration diﬀerences for each Nglobal Nglobal ∈ {5, 6, 7, 8, 9, 10}, at each numeric-volume, vvolume (xj ), for grid points x1 and x4 and for grid points x2 and x3 . Figure 2.1 also shows log plots of numerical concentration diﬀerences for each Nglobal ∈ {5, 6, 7, 8, 9}, at numeric-volume diﬀer10 n ences, |vvolume (xj ) − vvolume (xj )|, n ∈ {5, 6, 7, 8, 9}, for grid points x1 and x4 and for grid points x2 and x3 . A linear fit to the log plot of concentration diﬀerences at volume diﬀerences shows the relationship: 10 n 10 n |fi,t (xj ) − fi,t (xj )| = 17.9 · |vvolume (xj ) − vvolume (xj )|1.1 , (2.161) at grid points x1 and x4 , and 10 n 10 n |fi,t (xj ) − fi,t (xj )| = 7.7 · |vvolume (xj ) − vvolume (xj )|0.9 , (2.162) at grid points x2 and x3 , n ∈ {5, 6, 7, 8, 9}. Therefore, concentration diﬀerences converge locally in numeric-volume diﬀerences at approximately order 1.0. Numericvolume diﬀerences are approximately the diﬀerence of cubed grid spacing, hence concentration diﬀerences converge locally in grid spacing at approximately order 3, above the minimum spatial convergence order of 2, as predicted by finite diﬀerence formula derivations. 40 j log |f10(x )−fn (x )| i,t 0.008 0.006 10 i,t i,t j i,t −2.2 −2.4 j j |f10(x )−fn (x )| (µm−3) −2 0.01 0.004 0.002 −2.6 −2.8 −3 −3.2 0 2 4 6 8 n vvolume(xj) 10 3 (µm ) 12 −4.2 −4 −4 −3.8 (a) −3.4 −3.2 −3 (b) −2 i,t 0.004 0.002 0 2 4 6 8 n vvolume(xj) 10 3 (µm ) i,t 0.006 −2.2 j 0.008 j j 0.01 −2.4 10 i,t j log |f10(x )−fn (x )| 0.012 i,t |f10(x )−fn (x )| (µm−3) −3.6 10 n log10|vvolume(xj)−vvolume(xj)| x 10 −2.6 −2.8 −3 −4.2 12 −4 −4 −3.8 −3.6 −3.4 −3.2 −3 10 n log10|vvolume(xj)−vvolume(xj)| x 10 (c) (d) Figure 2.1: Spatial convergence of the base 2 multi-time-step forward Euler method in rod-cell geometry. Numerical concentration diﬀerences, for each Nglobal ∈ {5, 6, 7, 8, 9}, at each numeric-volume, are shown for grid points x1 and x4 in (a) and for grid points x2 and x3 in (c). Logs of numerical concentration diﬀerences, for each Nglobal ∈ {5, 6, 7, 8, 9}, at logs of grid point volume diﬀerences, and linear fit, are shown for grid points x1 and x4 in (b) and for grid points x2 and x3 in (d). n ∈ {5, 6, 7, 8, 9}. 41 2.4.7 Temporal Convergence of the Multi-Time-Step Forward Euler Method Numerical solutions to the system of diﬀusion equations, Equation 2.2, in rod-cell geometry, using the base 2 multi-time-step forward Euler method, should converge temporally to an exact solution, as time steps decrease in size to zero. An ordered vector containing numerical solutions, of reaction species, i, at time t, under macro ∆t time step, ∆t, and time step fine-tuning parameter, αi = 1, is denoted as fi,t . To ∆t test for temporal convergence, numerical solutions, fi,t , were found over the time interval, 10−3 s, using macro time steps, ∆tn = 2−n ·10−3 s, for all n ∈ {0, 1, 2, 3, 4, 5}, when starting from initial condition ftest , shown in Section 2.4.6, under a grid with Nglobal = 5, in a rod-cell of body length, L = .5µm, radius, R = .5µm, and diﬀusion coeﬃcient, Di = 1µm2 s−1 . Diﬀerences in numerical solutions from the numerical ∆t5 ∆tn solution found under macro time step, ∆t5 , ||fi,t − fi,t ||, were calculated for each macro time steps, ∆tn , n ∈ {0, 1, 2, 3, 4, 5}. Results are shown in Figure 2.2. A linear fit to a log of diﬀerences in numerical solutions from the numerical solution found under macro time step, ∆t5 , at the log of diﬀerences in macro time step from the macro time step, ∆t5 , shown in Figure 2.2, shows the relationship: ∆t5 ∆tn ||fi,t − fi,t || = 2.5 · 103 · |∆t5 − ∆tn |1.0 , (2.163) for all n ∈ {0, 1, 2, 3, 4}. Hence, numerical solutions converge temporally at approximately O(∆t), as expected for a forward Euler method. 42 || −3 −f∆i,t t=2 −3 −2 log2||f∆i,t t=2 ||f∆i,t t=2 0.5 0 −15 −14 −13 −12 −11 −n −10 −9 −3 1 −1 −5 1 −5 ⋅10 −3 −f∆i,t t=2 1.5 2 ⋅10 −n ⋅10 || 2 −n ⋅10 −3 2.5 0 −3 −4 −5 −15 −14 −13 −5 log2(∆ t=2 ⋅10 ) −12 −3 −11 −5 −10 −n log2|∆ t=2 ⋅10 −∆ t=2 ⋅10 | (a) (b) Figure 2.2: Temporal convergence of the base 2 multi-time-step forward Euler method. Diﬀerences in numerical solutions found under all macro time steps, ∆tn = 2−n · 10−3 s, n ∈ {0, 1, 2, 3, 4, 5}, from the numerical solution found under macro time step, ∆t5 , at the log of each macro time step, are shown in (a). Logs of the diﬀerences in numerical solutions found under all macro time steps, ∆tn , n ∈ {0, 1, 2, 3, 4}, from the numerical solution found under macro time step, ∆t5 , at the log of diﬀerences in macro time step from the macro time step, ∆tn , and linear fit are shown in (b) 2.4.8 Convergence of the Multi-Time-Step Forward Euler Method to the Infinite domain solution The known solution to the diﬀusion equation in three dimensions, on an infinite domain, starting with δ function at (x = x0 , y = y0 , z = z0 , t = 0) is: � � 1 (x − x0 )2 + (y − y0 )2 + (z − z0 )2 u(x, y, z, t) = exp − . (2.164) (4πDt)3/2 4Dt Numerical convergence studies demonstrate both spatial, in Section 2.4.6, and temporal, in Section 2.4.7, convergence of numerical solutions to the system of diﬀusion equations, Equation 2.2, in rod-cell geometry, using the base 2 multi-time-step forward Euler method. Beyond spatial and temporal convergence, numerical diﬀusion solutions must converge to the actual solution, which is approximately u(x, y, z, t), equation 2.164, when starting from a gaussian initial condition, centered at a grid 43 point suﬃciently far from the rod-cell membrane, with suﬃciently small variance, such that numerical solutions at rod-cell boundary grid points have suﬃciently small values, for some suﬃciently small time interval. In this regime, numerical solutions on the bounded rod-cell domain should closely approximate the infinite-domain solution. To test convergence of numerical solutions to u(x, y, z, t), a numerical solution was found using the suﬃciently small macro time step, ∆t = 10−4 s, over the suﬃciently small time interval, [0, 5 · 10−3 s], starting from the initial condition u(x, y, z, t0 = 4 · 10−3 s), centered at the interface of the rod-cell bottom cap and body, suﬃciently far from the the rod-cell boundary, at the cartesian coordinate representation of the grid point with r-index = 1, θ-index = 1, and z-index = 1, under a grid with Nglobal = 10, time step fine-tuning parameter, αi = 1, in a rod-cell of body length, L = .5µm, with radius, R = .5µm, and diﬀusion coeﬃcient, D = Di = 1µm2 s−1 . If the numerical solution is a good approximation to u(x, y, z, t), then the numerical solution should have a variance about the center of the initial Gaussian that grows like D · t. Hence, the variance, σ(t), was fit by nonlinear least squares to the numerical solution at each time, n · ∆t, n ∈ {0, 1, 2, . . . , 50}, using the MATLAB function “nlinfit”, shown in Figure 2.3. The resulting relationship is linear and well approximated by: σ(t) = 1.0 · t + 4.0 · 10−3 , (2.165) which is consistent with the expected linear spread in u, of rate 1s−1 , and initial time, t0 = 4 · 10−3 s. Relative diﬀerences between the numeric solution at each time, n · ∆t, n ∈ {0, 1, 2, . . . , 50}, and u(x, y, z, t0 + n · ∆t) should be small if numerical solutions are converging accurately. Denoting ut0 +n·∆t as an ordered vector, of which each entry is a function value, u(x, y, z, t0 + n · ∆t)), at a specific grid point, the global relative error, at each time, is then calculated: ||ut0 +n·∆t − fi,n·∆t || , ||ut0 +n·∆t || (2.166) and was found to monotonically increase in time, attaining a maximal value less than 6 · 10−3 , shown in Figure 2.3. A similar spread in the base 2 multi-time-step forward Euler diﬀusion equation solution, in rod-cell geometry, to the infinite domain solutions, u, with small global error, provides numerical support for convergence to the the infinite domain solution, under suﬃciently consistent conditions. Additionally, similar spreads with similar relative errors occur when initial distributions are 44 centered at any point of r index 1, in both caps and throughout the body (results not shown). −3 −3 x 10 3 ||⋅||u 7 σ (t) 5 0 7.5 t +n⋅∆ t 8 i,n⋅∆ t 8.5 6.5 −f 6 t +n⋅∆ t 5.5 4 2 1 ||u 0 5 4.5 0 x 10 ||−1 6 1 2 3 n⋅∆ t (s) 4 0 0 5 1 2 3 n⋅∆ t (s) −3 x 10 (a) 4 5 −3 x 10 (b) Figure 2.3: Nonlinear fit the numerical solutions to the infinite domain solution. (a) shows the nonlinear least squares fit of function u(x, y, z, t), in unknown timedependent parameter, σ(t), to the numerical solution at each time, n · ∆t. (b) shows the global relative error at each time, n · ∆t. n ∈ {0, 1, 2, . . . , 50}, ∆t = 10−4 s. 2.4.9 Diﬀusion on the Rod-Cell Boundary Biological systems in rod-cell geometry, such as the Min system, employ membrane bound components. Models of systems in rod-cell geometry with non-zero membrane bound component diﬀusion require numerical diﬀusion solutions that are restricted to the rod-cell boundary. Within the restriction of the system of diﬀusion equations, Eq. 2.2, to the rod cell boundary, there is no diﬀusion in the r direction, and diﬀusion in other spatial variables only occurs at the rod cell boundary, where r = R. Surface diﬀusion is described by: � � ∂ f¯i (R, θ, z, t) 1 ∂ 2 f¯i (R, θ, z, t) ∂ 2 f¯i (R, θ, z, t) = Di + , (2.167) ∂t R2 ∂θ2 ∂z 2 on the rod-cell body boundary, � � ∂¯ gi (R, θ, φ, t) 1 ∂ 2 g¯i (R, θ, φ, t) cos φ ∂¯ gi (R, θ, φ, t) 1 ∂ 2 g¯i (R, θ, φ, t) =Di + 2 + 2 , ∂t ∂θ2 R sin φ ∂φ R ∂φ2 R2 sin2 φ (2.168) 45 on the bottom cap boundary, and � ¯ i (R, θ, φ, t) ¯ i (R, θ, φ, t) ¯ i (R, θ, φ, t) ¯ i (R, θ, φ, t) � ∂h 1 ∂ 2h cos φ ∂ h 1 ∂ 2h =Di + 2 + 2 , ∂t ∂θ2 R sin φ ∂φ R ∂φ2 R2 sin2 φ (2.169) on the top cap boundary, for each membrane-bound reaction species i. Boundary conditions and cap-body interface matching conditions for diﬀusion on the rod-cell boundary are a restriction, to r = R, of boundary conditions for diﬀusion in rod-cell geometry. Numerical diﬀusion solutions on the rod-cell boundary are found using the rod-cell grid, as defined in Section 2.4, pertinent rod-cell finite diﬀerence formulas, as defined in Section 2.4.1, and pertinent rod-cell boundary and interface matching conditions, as defined in Section 2.4.2. Time-stepping follows the base 2 multi-timestep forward Euler method, as described in Section 2.4.3. Thus, base 2 multi-timestep forward Euler finite diﬀerences for diﬀusion on the rod-cell boundary: � � ¯ i (rj , θk , φm , tn + ∆τs ) f¯i (rj , θk , zl , tn + ∆τs ) + g¯i (rj , θk , φm , tn + ∆τs ) + h � � ¯ i (rj , θk , φm , tn ) = − f¯i (rj , θk , zl , tn ) + g¯i (rj , θk , φm , tn ) + h �� � � δl1 + δlNz 1 ∆τs · Di 1 − · δjNr · χ∆τs (Di · rj−2 ) · 2 ∆2θ f¯i (rj , θk , zl , tn ) 2 rj + (δl1 + δlNz ) · χ∆τs (Di · rj−2 ) · ∆2z f¯i (rj , θk , zl , tn ) �� 2¯ + (1 − δl1 − δlNz ) · χ∆τs (Di ) · ∆z fi (rj , θk , zl , tn ) �� � � δmNφ 1 +∆τs · Di 1− · δjNr · χ∆τs (Di · rj−2 sin−2 φm ) · 2 2 ∆2 g¯i (rj , θk , φm , tn ) 2 rj sin φm θ cos φm , rj−2 }) · 2 rj sin φm cos φm + χ∆τs (Di · max{ 2 , r−2 }) · rj sin φm j + χ∆τs (Di · max{ 46 cos φm ∆1φ g¯i (rj , θk , φm , tn ) 2 rj sin φm �� 1 2 ∆ g¯i (rj , θk , φm , tn ) rj2 φ +∆τs · Di �� δmNφ 1− 2 � · δjNr · � χ∆τs (Di · rj−2 sin−2 φm ) · cos φm , rj−2 }) · 2 rj sin φm cos φm , r−2 }) · + χ∆τs (Di · max{ 2 rj sin φm j + χ∆τs (Di · max{ rj2 1 ¯ i (rj , θk , φm , tn ) ∆2θ h 2 sin φm cos φm ¯ i (rj , θk , φm , tn ) ∆1φ h 2 rj sin φm �� 1 2¯ ∆ hi (rj , θk , φm , tn ) rj2 φ (2.170) F¯i,s is the Ngrid × Ngrid sparse matrix, such that: F¯i,s¯fi,t = ∆s¯fi,t , (2.171) where ¯fi,t , is a column vector with Ngrid entries that correspond to the concentration of membrane bound reaction species i, at time t, at each grid point, ordered from 1 to Ngrid , and ∆s¯fi,t is a column vector with Ngrid entries that correspond to the rod-cell boundary base 2 multi-time-step forward Euler finite diﬀerence value: � � ¯ i (rj , θk , φm , tn + ∆τs ) f¯i (rj , θk , zl , tn + ∆τs ) + g¯i (rj , θk , φm , tn + ∆τs ) + h � � ¯ i (rj , θk , φm , tn ) , − f¯i (rj , θk , zl , tn ) + g¯i (rj , θk , φm , tn ) + h for the time step, ∆τs , of membrane-bound reaction species, i, at time, t, at each grid point, ordered from 1 to Ngrid . A¯i is the complete rod-cell boundary base 2 multi-time-step forward Euler matrix, covering all finite diﬀerence evolutions over the macro time step, ∆t. A¯i satisfies: ¯fi,t+∆t = A¯i¯fi,t , (2.172) and is calculated: � A¯i = . . . �� � I + F¯i,1 �2 + F¯i,2 �2 + F¯i,3 �2 ... �2 2 + F¯i,Nt−1 + F¯i,Nt , (2.173) where I is a Ngrid × Ngrid sparse identity matrix. If A˘i = A¯i , excluding all dimensions that do not correspond to grid points on the rod-cell boundary, then as shown for rod-cell diﬀusion in Section 2.4.5, conservation during each time step of numerical diﬀusion, on the boundary, is determined by eigenvalue, 1, with corresponding strictly positive or strictly negative eigenvector, 47 ˘ , of A˘Ti . If 1 is a unique eigenvalue of A˘Ti , then an ordered vector, of which each v element is the numeric-surface area of a single grid-point on the rod-cell boundary, vsurface area , is calculated: vsurface area = rod-cell surface area ˘, v �Nθ (2Nφ +Nz −2) v p p=1 (2.174) ˘. where v˘p is the pth element of v The MATLAB function “eigs” was used to verify unique eigenvalue, 1, and corre˘ , of A˘Ti . sponding strictly positive or strictly negative eigenvector, v Numerical solutions to the system of diﬀusion equations, Equation 2.2, in rodcell geometry, using the base 2 multi-time-step forward Euler method, have shown convergence spatially, temporally, and to the infinite domain solution, u(x, y, z, t), Equation 2.164, under conditions stated in Section 2.4.8. Diﬀusion on the rod-cell boundary is a restriction of rod-cell diﬀusion, lacking only diﬀusion in the r direction, hence should be converging in all spatial dimensions, θ, z, and φ, thus spatial convergence, temporal convergence, and convergence to infinite domain solution are not shown, but assumed. 2.5 Reactions in Rod-Cell Geometry Reaction components, F (f (x, t)), in the system of reaction diﬀusion equations, Equation 2.1, are numerically solved at each grid point in rod-cell geometry using the Runge-Kutta 4 method, over the macro time step, ∆t. At each macro time step, the change from each reaction is added to numerical diﬀusion solutions to form a complete solution to the reaction-diﬀusion system, Equation 2.1. This step splitting procedure is reasonable if reactions are relatively slow compared to the diﬀusive time scale. Thus, macro time steps are chosen to be suﬃciently short, such that concentration changes from reactions are not significantly greater than concentration changes from diﬀusion, over the macro time step. 48 Chapter 3 Numerical Solutions to Reaction-Diﬀusion Equations on Changing Domains in Rod-Cell Geometry In Eschericia Coli, the Min system is tightly coupled with cell growth and division. Numerical methods are described to solve reaction diﬀusion equations in growing and dividing rod-cell geometries. 3.1 Previously Published Results and Numerical Methods for Solving Reaction-Diﬀusion Systems in Changing Domains Changing domains may drive specific biological pattern formation. On the body of the angelfish Pomacanthus, new stripes emerge between older, segregating stripes, during growth, a phenomena explained by Kondo and Asai’s reaction-diﬀusion system in a one dimensional growing domain (Kondo and Asai, 1995), and by Painter et al. (Painter et al., 1999) and Varea et al. (Varea et al., 1997) in a two dimensional growing domain. Kulesa et al. demonstrated the spatial patterning of alligator Teeth Primordia, through a one dimensional reaction-diﬀusion mechanism, in an exponentially growing domain, representing jaw growth (Kulesa et al., 1996). Crampin et al. explored a non-specific Turing type system in one spatial dimension 49 using various types of domain growth, and showed that their system may develop indefinite solution mode doubling under exponential domain growth, mode doubling with terminally sustained modes under logistic domain growth, and mode doubling that develops then“breaks up” under linear domain growth (Crampin et al., 1999). Modeling the Min system in vivo, Meinhardt and de Boer demonstrated mode splitting in a growing one dimensional domain (Meinhardt and de Boer, 2001). Numerically, Varea et al. and Kulesa et al. expanded their domains by incrementally stretching the spacing between grid points. Crampin et al., Kondo and Asai, Painter et al., and Meinhardt and de Boer enlarged their domains by incrementally adding new grid points with interpolated solution values. In this chapter, I develop methods for three dimensional, axial, rod cell growth by adding grid points, stretching the spacing between grid points, and simultaneously adding grid points and stretching the spacing between grid points. Several groups have modeled Min system behavior in dividing Eschericia Coli cells. Ventura and Sourjik solved their stochastic version of the Huang et al. model (Huang et al., 2003), under various static septum widths to postulate a possible mechanism for Min protein distribution in daughter cells (Ventura and Sourjik, 2011). Tostevin and Howard’s stochastic, one dimensional, Min system model accomplished dynamic septation by sequentially reducing the probability of diﬀusion across the domain mid point (Tostevin and Howard, 2006). Sengupta and Rutenberg deterministically modeled dynamic partitioning between daughter cells, using the Huang et al model, in a three dimensional cylindrical representation of Eschericia Coli, through sequential closure of diﬀusion across the axial midpoint, at radial grid points, with independent membrane interactions on each side of the growing partition (Sengupta and Rutenburg, 2007). In this chapter, I develop a numerical method for rod-cell division, similar to the Sengupta and Rutenberg method, though permitting diﬀusion of membrane-bound reaction species, with variable septation locations. 3.2 Numerically Solving Reaction-Diﬀusion Equations in Growing Rod-Cell Geometry Wild type Eschericia Coli cells grow axially, by adding peptidoglycan to the cell wall (Huang et al., 2008). Numerically, rod-cells may grow axially, by incrementally adding cylindrical volumes within a specific location in the rod-cell body, or by 50 incrementally stretching the rod-cell body. 3.2.1 Rod-Cell Growth by Adding Cylindrical Sections For rod-cell growth by adding cylindrical sections within the rod-cell body, first, a growth location in z is specified, zg , at z-index g ∈ {2, . . . , Nz − 1}. Then, concentrations of each reaction species, i, at grid points with z-index greater than g, are shifted in z-index by 1. The rod-cell has now grown by ∆z, and contains a cylinder of length ∆z, centered at zg+1 , which is empty in each reaction species, i. To fill the empty cylinder, for each reaction species, i, the value of the reaction species, at each grid point in zg+1 , is linearly interpolated from the values of the reaction species at grid points of equal r and θ locations in zg and zg+2 . Then, reaction-diﬀusion solutions are found using the base 2 multi-time-step forward Euler and RK4 methods, as described in Chapter 2. After any number of macro time steps, n · ∆t, n ∈ {0, 1, 2, . . . }, the rod-cell may be grown again, using any zg , g ∈ {2, . . . , Nz − 1}. 3.2.2 Rod-Cell Growth by Stretching the Rod-Cell Body For rod-cell growth by body stretching, first, a rod-cell final body length, LN ∈ {.5, 1, 1.5, 2, . . . }, is set. Then, Nz = 2LN Nglobal + 1, is set, and will be used for all rod-cell body lengths during growth. Next, reaction-diﬀusion equations are numerically solved, using the base 2 multi-time-step forward Euler and RK4 methods, as described in Chapter 2, under the initial body length, L0 . After any number of macro time steps, n · ∆t, n ∈ {0, 1, 2, . . . }, the rod cell grows by length ∆L1 . Growth occurs by stretching ∆z0 , the value of grid point spacing in z under body length L0 , to ∆z1 , the value of grid point spacing in z under body length L1 : L0 , Nz − 1 L1 L0 + ∆L1 ∆z1 = = . Nz − 1 Nz − 1 ∆z0 = (3.1) (3.2) The value of each reaction species, i, is unchanged at every grid point, but the rodcell has grown in length. Reaction-diﬀusion equations are then numerically solved, using the base 2 multi-time-step forward Euler and RK4 methods, under the rodcell body length L1 . The process is repeated until the final rod-cell body length, LN , has been attained. Then, reaction-diﬀusion solutions may be calculated for the 51 rod-cell, of body length LN , for any wanted time interval. To ensure the final body length, LN , ∆Lp , p ∈ {1, . . . , N }, must satisfy: LN − L 0 = 3.2.3 N � ∆Lp (3.3) p=1 A Comparison of Rod-Cell Growing methods Each method of rod-cell growth has benefits and drawbacks. Growth by adding cylindrical sections within the rod-cell body is computationally more eﬃcient than growth by body stretching, as it employs fewer grid point calculations, but requires strict growth by ∆z at each growth step, while growth by body stretching occurs through any chosen growth length, ∆Lp . Thus, growth by body stretching may implement smaller growth lengths than growth by adding cylindrical sections, resulting in smoother growth. Each growth method spatially distorts solutions. Growth by adding cylindrical sections within the rod-cell body only spatially distorts solutions where a growth cylinder has been added, but spatial distortion across the growth cylinder is relatively large, of O(∆z). Growth by body stretching spatially distorts solutions at all grid points in the rod-cell body, but spatial distortion at each grid point is p relatively small, of O( N∆L ). z −1 3.2.4 Hybrid Growth, by Adding Cylindrical Sections and Stretching Hybrid growth, growth by both adding cylindrical sections to the rod-cell body and stretching the rod-cell body, is as computational eﬃcient as growth by adding cylindrical sections to the rod-cell body, while permitting growth by any chosen growth length during each growth step, for smoother growth, as in growth by stretching the rod-cell body. Over each indexed macro growth step, q, which spans rod-cell growth by ∆z, a single cylindrical section is added to the rod-cell body; all other growth during q occurs by rod-cell body stretching. At each growth step p, p ∈ {1, . . . , Nq }, within each macro growth step q, the rod cell grows by ∆Lqp . During hybrid growth, Nz is set, as during growth by adding cylindrical sections to the rod-cell body. A growth amount, ∆Lq1 , is chosen. Then, to preserve equal spacing in z for growth by adding a cylindrical section, the rod-cell is stretched by .5∆Lq1 . A cylindrical section, of length .5∆Lq1 , is added to the rod-cell body at growth location, zg , 52 g ∈ {2, . . . , Nz − 1}. Thus, the rod-cell has grown by ∆Lq1 . Reaction-diﬀusion equations are then numerically solved, using the base 2 multi-time-step forward Euler and RK4 methods, as described in Chapter 2. After any number of macro time steps, n · ∆t, n ∈ {0, 1, 2, . . . }, the rod cell grows agains. For each consecutive growth step, p ∈ {2, . . . , Nq }, the rod-cell body is stretched by ∆Lqp . Between each growth step, p, reaction diﬀusion equations are solved for any chosen number of macro time steps. After the final growth step, Nq , under the macro growth step, q, the rod cell has grown by ∆z, hence ∆Lqp must be chosen such that: ∆z = Nq � ∆Lqp . (3.4) p=1 The process then continues for each successive macro growth step, q, until reaching the desired final cell length. During hybrid growth, by both adding cylindrical sections to the rod-cell body and stretching the rod-cell body, spatial distortion of solutions, at each grid ∆Lqp point within the rod-cell body, is of O( Nz −1 ), . 3.3 Dividing Rod-Cells During rod-cell division, an infinitely thin disc is introduces perpendicular to the z axis, within the rod cell body, which acts as a numeric invagination of the rod-cell boundary that successively closes in time, ultimately resulting in two independent rod-cell sub-domains. During the rod-cell division process, a division site in z, zd , d ∈ {2, . . . , Nz − 1}, is chosen. An infinitely thin disc, of radius, R, centered at r = 0, is inserted at zd . Each solution value, at each grid point of z-index greater than d, is shifted in z-index by 1. Shifted solution values maintain the same z location throughout the shift in z-index. The value of each solution, at each grid point within zd , is duplicated at the grid point of equal r and θ location in zd+1 . Thus, an infinitely thin disk has been inserted at zd , with equal solution values on opposite sides, at zd and zd+1 , which will act as boundary points, separating the two rod-cell sub-domains as the boundary invaginates. The rod-cell boundary successively invaginates equally along both sides of the infinitely thin disk, in the r direction, from r = R to r = 0, acting as a rod-cell diﬀusion barrier for bulk reaction species and as a domain extension for membrane 53 bound reaction species. Initially, there is no invagination. Then, at each invagination step, q ∈ {1, 2, . . . , Nr }, the boundary invaginates to R − (2(q−1)+1)∆r , the 2 distance half-way between grid points of r-indices Nr − (q − 1) and Nr − q, for q ∈ {1, 2, . . . , Nr − 1}, and r = 0, for q = Nr . To ensure smooth transitions between invagination steps, the value of each membrane-bound reaction species, at each newly invaginated grid point, takes on the value of the previously invaginated grid point, of equal θ orientation and z-index. Before the first invagination step, between any successive invagination steps, and after the final invagination step, any number of macro time steps, n · ∆t, n ∈ {0, 1, 2, . . . } may be taken. 3.3.1 The Base 2 Multi-Time-Step Forward Euler Method for Solving Diﬀusion Equations in Dividing Rod-Cells Solutions to the diﬀusion parts of reaction-diﬀusion equations, in dividing rod-cells, are found using the base 2 multi-time-step forward Euler method, described in Chapter 2. All spatial finite diﬀerence formulas, from Section 2.4.1, and all boundary conditions, from Section 2.4.2, for solving diﬀusion equations in rod-cell geometry are applicable for solving diﬀusion equations in dividing rod-cells, but extra conditions are required. Grid points at zd are separated from grid points at zd+1 by an infinitely thin disc. Hence, grid points within zd have a measure 0 diﬀerence in location from grid points within zd+1 , when equally located in the r and θ directions. Thus, until grid points are separated by the invaginated boundary, solution values at grid points within zd should remain identical to solution values within zd+1 , at grid points of equal r and θ locations. To ensure equal solution values at grid points of equal r and θ locations, at zd and zd+1 , when diﬀusion is open in z across zd and zd+1 , while maintaining the proper magnitude of diﬀusivity, the total diﬀusion in r and θ, at zd and zd+1 , is half the diﬀusion in r and θ at zd , and half the diﬀusion in r and θ at zd+1 . At each invagination step, q, the boundary has invaginated in the r direction past grid points of r-index Nr − (q − 1). Thus, at grid points of r-index ≥Nr − (q − 1), diﬀusion is closed in the positive z direction at zd , and closed in the negative z-direction at zd+1 . 54 Second order finite diﬀerences in z, at zd and zd+1 are: ∂ 2 f (rj , θk , zd , tn ) f (rj , θk , zd−1 , t) 2f (rj , θk , zd , t) f (rj , θk , zd+2 , t) ≈ − + , ∂z 2 ∆z 2 ∆z 2 ∆z 2 (3.5) ∂ 2 f (rj , θk , zd+1 , tn ) f (rj , θk , zd−1 , t) 2f (rj , θk , zd+1 , t) f (rj , θk , zd+2 , t) ≈ − + . (3.6) ∂z 2 ∆z 2 ∆z 2 ∆z 2 To impose closed boundaries in the positive z direction across zd : fi (rj , θk , zd+2 ) = fi (rj , θk , zd ), (3.7) in Eq. 3.5, for j ≥ Nr − (q − 1). To impose closed boundaries in the negative z direction actross zd+1 : fi (rj , θk , zd−1 ) = fi (rj , θk , zd+1 ), (3.8) in Eq. 3.6, for j ≥ Nr − (q − 1). Additional conditions required for solving diﬀusion equations in dividing rodcells, in conjunction with the base 2 multi-time-step forward Euler method for solving diﬀusion equations in rod-cell geometry, from Section 2.4.3, result in base 2 multi-time-step forward Euler finite diﬀerences for solving diﬀusion equations in dividing rod-cells: (fi (rj , θk , zl , tn + ∆τs ) + gi (rj , θk , φm , tn + ∆τs ) + hi (rj , θk , φm , tn + ∆τs )) − (fi (rj , θk , zl , tn ) + gi (rj , θk , φm , tn ) + hi (rj , θk , φm , tn )) = �� � δl1 − δlNz ∆τs · Di 1 − − δld − δl(d+1) · 2 � χ∆τs (Di · max{1, rj−1 }) · ∆2r fi (rj , θk , zl , tn ) + χ∆τs (Di · max{1, rj−1 }) · + χ∆τs (Di · rj−2 ) · 1 1 ∆ fi (rj , θk , zl , tn ) rj r 1 2 ∆ fi (rj , θk , zl , tn ) rj2 θ + (δl1 + δlNz ) · χ∆τs (Di · rj−2 ) · ∆2z fi (rj , θk , zl , tn ) � + (1 − δl1 − δlNz ) · χ∆τs (Di ) · ∆2z fi (rj , θk , zl , tn ) 55 � δld + · (1 + H[j − (Nr − (q − 1))]) · χ∆τs (Di · max{1, rj−1 }) · ∆2r fi (rj , θk , zl , tn ) 2 1 + χ∆τs (Di · max{1, rj−1 }) · ∆1r fi (rj , θk , zl , tn ) rj � 1 2 −2 + χ∆τs (Di · rj ) · 2 ∆θ fi (rj , θk , zl , tn ) rj � δld + · (1 − H[j − (Nr − (q − 1))]) · χ∆τs (Di · max{1, rj−1 }) · ∆2r fi (rj , θk , zl+1 , tn ) 2 1 + χ∆τs (Di · max{1, rj−1 }) · ∆1r fi (rj , θk , zl+1 , tn ) rj � 1 2 −2 + χ∆τs (Di · rj ) · 2 ∆θ fi (rj , θk , zl+1 , tn ) rj � δl(d+1) + · (1 + H[j − (Nr − (q − 1))]) · χ∆τs (Di · max{1, rj−1 }) · ∆2r fi (rj , θk , zl , tn ) 2 1 + χ∆τs (Di · max{1, rj−1 }) · ∆1r fi (rj , θk , zl , tn ) rj � 1 2 −2 + χ∆τs (Di · rj ) · 2 ∆θ fi (rj , θk , zl , tn ) rj � δl(d+1) + · (1 − H[j − (Nr − (q − 1))]) · χ∆τs (Di · max{1, rj−1 }) · ∆2r fi (rj , θk , zl−1 , tn ) 2 1 + χ∆τs (Di · max{1, rj−1 }) · ∆1r fi (rj , θk , zl−1 , tn ) rj � 1 2 −2 + χ∆τs (Di · rj ) · 2 ∆θ fi (rj , θk , zl−1 , tn ) rj � � fi (rj , θk , zl−1 , tn ) 2fi (rj , θk , zl , tn ) fi (rj , θk , zl+2 , tn ) + δld · χ∆τs (Di ) · − + ∆z 2 ∆z 2 ∆z 2 � �� fi (rj , θk , zl−2 , tn ) 2fi (rj , θk , zl , tn ) fi (rj , θk , zl+1 , tn ) + δl(d+1) · χ∆τs (Di ) · − + ∆z 2 ∆z 2 ∆z 2 56 +∆τs · Di �� δmNφ 1− 2 � � · χ∆τs (Di · max{1, 2rj−1 }) · ∆2r gi (rj , θk , φm , tn ) 2 1 ∆ gi (rj , θk , φm , tn ) rj r 1 + χ∆τs (Di · rj−2 sin−2 φm ) · 2 2 ∆2θ gi (rj , θk , φm , tn ) rj sin φm cos φm cos φm + χ∆τs (Di · max{ 2 , rj−2 }) · 2 ∆1 gi (rj , θk , φm , tn ) rj sin φm rj sin φm φ �� 1 2 cos φm −2 , r }) · 2 ∆φ gi (rj , θk , φm , tn ) + χ∆τs (Di · max{ 2 rj sin φm j rj �� � � δmNφ +∆τs · Di 1 − · χ∆τs (Di · max{1, 2rj−1 }) · ∆2r hi (rj , θk , φm , tn ) 2 2 + χ∆τs (Di · max{1, 2rj−1 }) · ∆1r hi (rj , θk , φm , tn ) rj 1 + χ∆τs (Di · rj−2 sin−2 φm ) · 2 2 ∆2θ hi (rj , θk , φm , tn ) rj sin φm cos φm cos φm + χ∆τs (Di · max{ 2 , rj−2 }) · 2 ∆1 hi (rj , θk , φm , tn ) rj sin φm rj sin φm φ �� cos φm 1 2 −2 + χ∆τs (Di · max{ 2 , r }) · 2 ∆φ hi (rj , θk , φm , tn ) rj sin φm j rj + χ∆τs (Di · max{1, 2rj−1 }) · (3.9) where H is the discrete Heaviside function. Using dividing rod-cell base 2 multi-time-step forward Euler finite diﬀerences, Ai , the complete base 2 multi-time-step forward Euler matrix, covering all finite diﬀerence evolutions, over the macro time step, ∆t, is calculated as in Section 2.4.3. Conservation during each time step of numerical diﬀusion, in the dividing rod-cell, is determined using the method of Section 2.4.5. Complete solutions to reactiondiﬀusion equation, in the dividing rod-cell, are found using the method of Section 2.5. 57 3.3.2 The Base 2 Multi-Time-Step Forward Euler Method for Solving Diﬀusion Equations on Dividing Rod-Cell Boundaries The dividing rod-cell boundary is essentially the same as the the rod-cell boundary, with added invagination in the r direction at zd and zd+1 . Solving diﬀusion equations on the dividing rod-cell boundary follows directly from solving diﬀusion equations on the rod-cell boundary, as in Section 2.4.9, with added conditions at invaginated boundary points. Membrane bound proteins can diﬀuse around the leading edge of the invaginated membrane, moving from one face of the invagination to the other, until the final invagination step, q = Nr . For each invagination step, q ∈ {1, 2, . . . , Nr−1 }, the boundary is chosen to invaginate half-way between grid points of r-indices, Nr − (q − 1) and Nr − q. Thus, grid points of r-index Nr − (q − 1), within zd and zd+1 , are separated in invaginated boundary by a distance, ∆r, when grid points are of equal θ orientation. Grid point locations in r, at equal θ orientation, relative to the leading edge of the invaginated boundary are shown in Figure 3.1. leading edge of invaginated boundary rN − (q−1) rN − (q−2) rN − (q−1) r r rN − (q−2) r r zd+1 zd Figure 3.1: Grid point locations in r relative to the leading edge of the invaginated boundary. Grid points, of equal θ orientation, on opposite faces of the invaginated boundary are shown. Separating grid points across the invaginated boundary by ∆r allows the implementation of symmetric five point finite diﬀerence formulas in r, shown in Section 58 2.4.1. To allow diﬀusion around the leading edge of the invaginated boundary, conditions: f¯i (rNr −q , θk , zd ) = f¯i (rNr −(q−1) , θk , zd+1 ), f¯i (rNr −q−1 , θk , zd ) = f¯i (rN −(q−2) , θk , zd+1 ), r f¯i (rNr −q , θk , zd+1 ) = f¯i (rNr −(q−1) , θk , zd ), f¯i (rNr −q−1 , θk , zd+1 ) = f¯i (rN −(q−2) , θk , zd ), r (3.10) (3.11) (3.12) (3.13) are set in finite diﬀerence formulas in r, for q ∈ {1, 2, . . . , Nr−1 }. When q = Nr , diffusion is closed across the invaginated boundary, hence finite diﬀerence conditions, in r, follow those for bulk rod-cell diﬀusion, shown in Section 2.4.2. As in the bulk of the dividing rod-cell, the dividing rod-cell boundary is closed in the positive z direction, for rj = R, across zd , and closed in the negative z direction across zd + 1, hence second order finite diﬀerence formulas, in z, and boundary conditions, at zd and zd + 1, follow those from Section 3.3.1, Equations 3.5, 3.6, 3.7, 3.8. Adding conditions for the invaginated boundary to base 2 multi-time-step forward Euler finite diﬀerences for solving diﬀusion equations on the rod-cell boundary, shown in 2.4.9, result in base 2 multi-time-step forward Euler finite diﬀerences for solving diﬀusion equations on the dividing rod cell boundary: � � ¯ i (rj , θk , φm , tn + ∆τs ) f¯i (rj , θk , zl , tn + ∆τs ) + g¯i (rj , θk , φm , tn + ∆τs ) + h � � ¯ i (rj , θk , φm , tn ) = − f¯i (rj , θk , zl , tn ) + g¯i (rj , θk , φm , tn ) + h �� � � 1 δl1 + δlNz ∆τs · Di 1 − − δld − δl(d+1) · δjNr · χ∆τs (Di · rj−2 ) · 2 ∆2θ f¯i (rj , θk , zl , tn ) 2 rj + (δl1 + δlNz ) · χ∆τs (Di · rj−2 ) · ∆2z f¯i (rj , θk , zl , tn ) � + (1 − δl1 − δlNz ) · χ∆τs (Di ) · ∆2z f¯i (rj , θk , zl , tn ) 59 � � + δld + δl(d+1) · H[j − (Nr − (q − 1))]· � χ∆τs (Di · max{1, rj−1 }) · ∆2r f¯i (rj , θk , zl , tn ) � 1 1¯ (3.14) + χ∆τs (Di · · ∆r fi (rj , θk , zl , tn ) rj �¯ fi (rj , θk , zl−1 , tn ) 2f¯i (rj , θk , zl , tn ) + δjNr · δld · χ∆τs (Di ) · − ∆z 2 ∆z 2 � f¯i (rj , θk , zl+2 , tn ) + ∆z 2 �¯ fi (rj , θk , zl−2 , tn ) 2f¯i (rj , θk , zl , tn ) + δjNr · δl(d+1) · χ∆τs (Di ) · − ∆z 2 ∆z 2 �� f¯i (rj , θk , zl+1 , tn ) + ∆z 2 �� � � δmNφ 1 +∆τs · Di 1− · δjNr · χ∆τs (Di · rj−2 sin−2 φm ) · 2 2 ∆2 g¯i (rj , θk , φm , tn ) 2 rj sin φm θ max{1, rj−1 }) cos φm cos φm , rj−2 }) · 2 ∆1 g¯i (rj , θk , φm , tn ) 2 rj sin φm rj sin φm φ �� cos φm 1 2 −2 + χ∆τs (Di · max{ 2 , r }) · 2 ∆φ g¯i (rj , θk , φm , tn ) rj sin φm j rj �� � � δmNφ 1 ¯ i (rj , θk , φm , tn ) +∆τs · Di 1− · δjNr · χ∆τs (Di · rj−2 sin−2 φm ) · 2 2 ∆2 h 2 rj sin φm θ + χ∆τs (Di · max{ cos φm , rj−2 }) · 2 rj sin φm cos φm + χ∆τs (Di · max{ 2 , r−2 }) · rj sin φm j + χ∆τs (Di · max{ cos φm ¯ i (rj , θk , φm , tn ) ∆1φ h 2 rj sin φm �� 1 2¯ ∆ hi (rj , θk , φm , tn ) rj2 φ (3.15) Using dividing rod-cell boundary base 2 multi-time-step forward Euler finite diﬀerences, A¯i , the complete base 2 multi-time-step forward Euler matrix, covering all finite diﬀerence evolutions, over the macro time step, ∆t, is calculated as in Section 2.4.9. Conservation during each time step of numerical diﬀusion, on the dividing rod-cell boundary, is determined using the method of Section 2.4.9. Complete solutions to reaction-diﬀusion equation, on the dividing rod-cell boundary, are found using the method of Section 2.5. 60 3.4 Corrections for Changes in Total Amounts of Conserved Quantities during Rod-Cell Growth and Division When rod-cells numerically grow or divide, the total amount of each conserved quantity is altered by changes in grid point numeric volumes or by interpolation at newly introduced interior and/or boundary grid points. Depending on the model, the total amount of each conserved quantity may be independent of cell length and volume, dependent on cell length, or dependent on cell volume. Concentrations may be added or subtracted to each reaction species, at each grid point, in any desired manner, such that the total amount of each conserved quantity, at each growth or division step, is corrected to be the desired total amount of each conserved quantity, based on the specific model. A specific example, using the Huang et al. model, is shown in Section 4.2. 61 Chapter 4 Examples of Reaction-Diﬀusion Solutions in Rod-Cell Geometry To demonstrate the methods, developed in Chapters 2 and 3, for numerically solving reaction-diﬀusion equations in rod-cell geometry, examples using the Huang et al. Min system model (Huang et al., 2003), in a representation of a wild type cell, which grows and divides, and in a representation of a filamentous mutant cell, which continuously grows, are shown. Any reaction-diﬀusion model may be solved separately under static cell size, growth, or division to develop a more in depth understanding of the model during each individual process. However, here, all conditions are shown together to demonstrate the more numerically complicated and biologically relevant solutions to the Huang et al. model. The examples shown are not intended to describe the Huang et al. model extensively, but to show the potential of the numerical methods described in previous chapters. Results and future applications are discussed. 4.1 The Huang et al. Model The Huang et al. model describes the in vivo behavior of Eschericia coli proteins, MinD and MinE. The model contains five states: ρD:ADP , the concentration of cytosolic MinD:ATP complexes, ρD:AT P , the concentration of cytosolic MinD:ATP complexes, ρE , the concentration of cytosolic MinE, ρd , the concentration of membrane bound MinD:ATP complexes, and ρde , the concentration of membrane-bound MinE:MinD:ATP complexes. In the model, cytosolic MinD:ADP complexes release 62 ADP →AT P ADP and bind an ATP, forming MinD:ATP complexes, at the rate, σD . Cytosolic MinD:ATP complexes bind to the membrane, at the background rate, σD , and are actively recruited to membrane by membrane-bound MinD:ATP and MinE:MinD:ATP complexes, at the rate, σdD . Once bound to the membrane, MinD:ATP complexes bind to cytosolic MinE, forming membrane-bound MinE:MinD:ATP complexes, at the rate, σE . Completing the reaction cycle, membrane bound MinE:MinD:ATP complexes undergo ATP hydrolysis, releasing MinD:ADP complexes and MinE into the cytosol, at the rate, σde . In the Huang et al. model, there is no diﬀusion on the boundary. Including diﬀusion on the boundary, the Huang et. al model becomes: ∂ρD:ADP ADP →AT P =DD ∇2 ρD:ADP − σD ρD:ADP + δ(r − R)σde ρde , ∂t ∂ρD:AT P ADP →AT P =DD ∇2 ρD:AT P + σD ρD:ADP ∂t − δ(r − R)[σD + σdD (ρd + ρde )]ρD:AT P , ∂ρE =DE ∇2 ρE + δ(r − R)σde ρde − δ(r − R)σE ρd ρE , ∂t ∂ρd =Dd ∇2 |R ρd − σE ρd ρE (R) + [σD + σdD (ρd + ρde )]ρD:AT P (R), ∂t ∂ρde =Dde ∇2 |R ρde − σde ρde + σE ρd ρE (R). ∂t (4.1) (4.2) (4.3) (4.4) (4.5) Parameters within the Huang et al. model have values: DD = DE =2.5µm2 s−1 , (4.6) Dd = Dde =0µm2 s−1 , (4.7) ADP →AT P σD −1 (4.8) σD =0.025µms−1 , (4.9) σdD =0.0015µm3 s−1 , (4.10) =1s , −1 σde =0.7s , (4.11) 3 −1 σE =0.093µm s . 63 (4.12) 4.2 The Huang et al. Model in a Wild Type Cell that Grows and Divides I numerically solved the Huang et al. model in a growing then dividing rod-cell, representing wild type Escherichia Coli behavior, under a grid with Nglobal = 5. The rod-cell grew by hybrid growth, adding cylindrical sections to the center of the rod-cell body and stretching the rod-cell body, as described in Section 3.2.4, from 1.5µm to 3µm in length, at a rate of 0.02µm every 12s. Once the rod-cell was 3µm in length, division began at the center of the rod-cell body, and continued at a uniform rate, over a time interval of 36s. Then, numerical solutions were calculated for an additional two minutes, to determine the behavior of the system in the divided rod-cell. Initial conditions are constant, chosen such that the total concentration of cytosolic MinD:ATP is 1000µm−1 and the total concentration of cytosolic MinE is 350µm−1 , with 25% added random noise at each grid point, sampled from a normal distribution, then scaled, using grid point volumes, to a total concentration of 1000µm−1 cytosolic MinD:ATP and 350µm−1 cytosolic MinE, as stipulated in the Huang et al. model. Solutions were calculated under the Huang et al reaction parameter, shown above. At each growth and division step, the total amount of MinD and MinE were corrected to maintain the same mean volumetric concentration as the initial condition. MinD total amounts were corrected by uniformly adding the necessary concentration of ρD:ADP , to each grid point, when positive MinD corrections were required, analogous to constant production of ρD:ADP . When negative MinD corrections were required, ρD:ADP and ρD:AT P were removed, at each grid point, proportional to the amount of each reaction species, analogous to constant degradation of cytosolic MinD proteins. MinE corrections followed the same method as MinD corrections. To visualize the solution, one-dimensional mean concentration projections were calculated, at each time step, by sorting grid points into bins based on axial location, then finding the mean concentration within each bin, using grid point volumes. Results are shown in Figure 4.1. 64 (a) (b) (c) (d) Figure 4.1: Huang et al. solutions in a wild type cell that grows and divides. Onedimensional mean solution concentrations are shown for each reaction species: (a), ρD:ADP , (b), ρD:AT P , (c), ρD:AT P , and (d), ρd . ρde resembles ρd in behavior, with a slight lag, hence is not shown. The rod-cell grows in length from 1.5µm to 3µm, over the time interval, [0s 913s). Rod-cell division occurs over the time interval, [913s 949s]. This example and variants (not shown) demonstrate several biological behaviors, providing insight into the Min system, and the Huang et al. model. Regular pole-to-pole oscillations do not begin until the rod-cell has attained a length of about 2µm, as shown in the wild type example. However, regular pole-to-pole oscillatory 65 behavior occurs at cell lengths less than 2µm when cytosolic diﬀusion coeﬃcients are less than 2.5µm2 s−1 , when starting from the same initial condition (results not shown), indicating that larger cytosolic protein diﬀusivity may connect reaction regions suﬃciently to disrupt regular pole-to-pole oscillations in short cells, a possible cause for the appearance of regular pole-to-pole oscillations only in cells longer than a specific length. In the wild type example, the minimum mean MinD concentration over time is located at the rod-cell center, while the maximum mean MinD concentrations over time are located at cell poles, as is observed experimentally, and has been demonstrated by models in static domains, (Loose et al., 2008), (Cytrynbaum and Marshall, 2007), (Huang et al., 2003), (Fischer-Friedrich et al., 2010), (Tostevin and Howard, 2006), (Kerr et al., 2006), (Fange and Elf, 2006), (Sengupta and Rutenburg, 2007), and in a one dimensional growing domain, (Meinhardt and de Boer, 2001). Regular pole-to-pole oscillations in MinD and MinE occur with periodicity of about 40s, within the range of quickly oscillating Min proteins, in vivo, (Huang et al., 2003). After rod-cell division, oscillations continue in one of rod-cell sub-domains, but not in the other, despite equal size. This is consistent with experimental observations, (Fischer-Friedrich et al., 2010), and model solutions, (Sengupta and Rutenburg, 2007), where one daughter cell contains regular pole-to-pole MinD oscillations and the other contains stochastic pole-to-pole MinD switching, despite similar sized daughter cells with similar MinD total amounts. This example shows MinE concentrations following MinD concentrations, during regular pole-to-pole oscillations. The lag in MinE behind MinD could result in an unequal MinE distribution amongst daughter cells, despite a roughly equal MinD distribution, resulting in diﬀering daughter cell behavior, a possible cause of the Fischer-Friedrich findings. No oscillatory behavior was observed when incorporating diﬀusion on the rod-cell boundary, when boundary diﬀusion coeﬃcients were chosen an order of magnitude smaller than cytosolic diﬀusion coeﬃcients (results not shown). Membrane bound proteins likely diﬀuse in vivo. A future investigation of how and to what extent diﬀusion of membrane bound proteins disrupt regular pole-to-pole oscillation in the Huang et al. model will provide a biological test of the model, and provide possible insight into model improvements. 66 4.3 The Huang et al. Model in a Continuously Growing Filamentous Mutant Cell I numerically solved the Huang et al. model in a continuously growing rod-cell, representing Escherichia Coli filamentous mutant behavior, under a grid with Nglobal = 5. The rod-cell was grown by hybrid growth, as described in Section 3.2.4, from 1.5µm to 10µm in length, at a rate of 0.02µm every 2s. Initial conditions, corrections for MinD and MinE total amounts, and one dimensional mean concentration projections were found as in the wild type example, Section 4.2. Results are shown in Figure 4.2 67 (a) (b) (c) (d) Figure 4.2: The Huang et al. model in a continuously growing filamentous mutant cell. One-dimensional mean solution concentrations are shown for each reaction species: (a), ρD:ADP , (b), ρD:AT P , (c), ρD:AT P , and (d), ρd . ρde resembles ρd in behavior, with a slight lag, hence is not shown. The rod-cell grows in length from 1.5µm to 10µm, over the time interval, [0s 853s]. During oscillations, in the filamentous mutant example, Min protein maximal concentrations (over time) are initially located only at the cell poles. As the cell increases in length, the Min protein pole-to-pole oscillation period increases. Then, a third maximal concentration (over time) appears near midcell, with Min protein pole-to-midcell oscillation periodicities similar to pole-to-pole oscillation periodici68 ties of the cell, when shorter. A similar result was demonstrated in a growing one dimensional domain by Meinhardt and de Boer (Meinhardt and de Boer, 2001). The formation of Min protein concentration maximums (over time), away from the poles, in long cells, resembles experimental findings from Raskin and de Boer (Raskin and de Boer, 1999), and model solutions in static domains, (Huang et al., 2003), (Tostevin and Howard, 2006), (Fange and Elf, 2006). A more complete comparison of Raskin and de Boer’s experimental findings with Huang et al. model solutions, in growing, three dimensional, filamentous cells, including the formation of multiple Min protein maximal sites over time, and changes is oscillation periodicity, require future investigation. 4.4 Future Applications Solving Min system models, using the computationally eﬃcient base 2 multi-timestep forward Euler method for solving reaction-diﬀusion equations in growing and dividing rod-cell geometries, will allow the study of models in terms of: the ability to develop regular pole-to-pole oscillatory solutions, as functions of the total amounts of each reaction species, diﬀusivity, and cell length; cell length dependent changes in solution periodicity, as rod-cells grow; the ability to undergo a bifurcation from a single pole-to-pole oscillation to two pole-to-midcell oscillations to higher order oscillations; the eﬀect of cell growth on bifurcation processes; and how divided rodcell sub-domain behaviors diﬀer, as functions of when the rod-cell divides, where the rod-cell divides, and how quickly the rod the-cell divides. Ultimately, a fairly complete description of various models’ behaviors, in the biologically relevant space of growing and dividing rod-cells, will develop a deeper understanding of mathematically interesting Min system models, provide insight into new models, and produce biologically pertinent predictions. 69 Bibliography E. J. Crampin, E. A. Gaﬀney, and P. K. Maini. Reaction and diﬀusion on growing domains: Scenarios for robust pattern formation. Bulletin of Mathematical Biology, 61:1093–1120, 1999. E. N. Cytrynbaum and D. L. Marshall. A multistranded polymer model explains MinDE dynamics in E. coli. Biophysical Journal, 93:1134–1150, aug 2007. P. A. de Boer, R. E. Crossley, and L. I. Rothfield. A division inhibitor and a topological specificity factor coded for by the minicell locus determine proper placement of the division septum in E. coli. Cell, 56(4):641–649, feb 1989. D. Fange and J. Elf. Noise-induced Min phenotypes in Escherichia coli. PLoS Computational Biology, 2:637–648, jun 2006. E. Fischer-Friedrich, G. Meacci, J. Lutkenhaus, H. Chat´e, and K. Karsten. Intraand intercellular fluctuations in Min protein dynamics decrease with cell length. PNAS, 107(14):6134–6139, apr 2010. K. C. Huang, Y. Meir, and N. S. Wingreen. Dynamic structures in Escherichia coli : Spontaneous formation of MinE rings in MinD polar zones. PNAS, 100(22): 12724–12728, oct 2003. K. C. Huang, R. Mukhopadhyay, B. Wen, and N. S. Wingreen. Cell shape and cell-wall organization in Gram-negative bacteria. PNAS, 105(49):19282–19287, 2008. V. Ivanov and K. Mizuuchi. Multiple modes of interconverting dynamic pattern formation by bacterial cell division proteins. PNAS, 107(18):8071–8078, may 2010. 70 R. A. Kerr, H. Levine, T. J. Sejnowski, and W.-J. Rappel. Division accuracy in a stochastic model of Min oscillations in Escherichia coli. PNAS, 103(2):347–352, jan 2006. S. Kondo and R. Asai. A reaction-diﬀusion wave on the skin of the marine angelfish Pomacanthus. Nature, 376:765–768, aug 1995. P. Kulesa, G. Cruywagen, S. Lubkin, P. Maini, J. Sneyd, M. Ferguson, and J. Murray. On a model mechanism for the spatial patterning of teeth primordia in the alligator. J. theor. Biol., 180:287–296, 1996. L. L. Lackner, D. M. Raskin, and P. A. J. de Boer. ATP-dependent interactions between Escherichia coli Min proteins and the phospholipid membrane in vivo. Journal of Bacteriology, pages 735–749, feb 2003. M.-C. Lai. A note on finite diﬀerence discretization for Poisson equation on a disk. Numer. Meth. Partial Diﬀ. Eq, 17:199–2003, 2001. M. Loose, E. Fischer-Friedrich, J. Ries, K. Kruse, and P. Schwille. Spatial regulators for bacterial cell division self-organize into surface waves in vitro. Science, 320: 789–792, may 2008. M. Loose, E. Fischer-Friedrich, C. Herold, K. Kruse, and P. Schwille. Min protein patterns emerge from rapid rebinding and membrane interaction of MinE. Nature Structural and Molecular Biology, 18(5):577–584, may 2011. G. Luyan Ma and L. F. King. Positioning of the MinE binding site on the MinD surface suggests a plausible mechanism for activation of the Escherichia coli MinD ATPase during division site selection. Molecular Microbiology, 54(1):99–108, oct 2004. H. Meinhardt and P. A. J. de Boer. Pattern formation in Escherichia coli : a model for the pole-to-pole oscillations of Min proteins and the localization of the division site. PNAS, 98(25):14202–14207, dec 2001. K. Painter, P. Maini, and H. Othmer. Stripe formation in juvenile Pomacanthus explained by generalized turing mechanism with chemotaxis. Proc. Natl. Acad. Sci. USA, 96:5549–5554, may 1999. 71 D. Raskin and P. de Boer. Rapid pole-to-pole oscillation of a protein required for directing division to the middle of Escherichia coli. Proc. Natl. Acad. Sci. USA, 96(9):4971–4976, apr 1999. S. Sengupta and A. Rutenburg. Modeling partitioning of Min proteins between daughter cells after septation in Escherichia coli. Phys. Biol., 4:145–153, 2007. T. H. Szeto, S. L. Rowland, C. L. Habrukowich, and G. F. King. The MinD membrane targeting sequence is a transplantable lipid-binding helix. J. Biol. Chem., 278:40050–40056, jul 2003. F. Tostevin and M. Howard. A stochastic model of Min oscillations in Escherichia coli and Min protein segregation during cell division. Physical Biology, 3:1–12, 2006. C. Varea, J. Arag´on, and R. Barrio. Confined turing patterns in growing systems. Physical Review E, 56(1):1250–1256, jul 1997. D. B. Ventura and V. Sourjik. Self-organized partitioning of dynamically localized proteins in bacterial cell division. Molecular Systems Biology, 7:457, 2011. H. Zonglin, E. P. Gogol, and J. Lutkenhaus. Dynamic assembly of MinD on phospholipid vesicles regulated by ATP and MinE. PNAS, 99(10):6761–6766, may 2002. 72
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A computationally efficient method for solving Min...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A computationally efficient method for solving Min system reaction-diffusion equations within growing… Carlquist, William Christopher 2012
pdf
Page Metadata
Item Metadata
Title | A computationally efficient method for solving Min system reaction-diffusion equations within growing and dividing domains that approximate a rod-shaped bacterium |
Creator |
Carlquist, William Christopher |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | Few biological systems, showing complex pattern formation that spans multiple spatial and temporal scales, have been reduced in understanding to several components. The Min system in Escherichia coli, consists of three proteins, MinD, MinE, and MinC. Through ordered, cyclic, membrane binding and unbinding, facilitated by ATP hydrolysis, the Min system regulates the site of cell division in vivo. The Min system is tightly coupled with cell growth and division. Various mathematical models have been proposed to describe specific biological phenomena, arising from the Min system, but no model has been tested in a biologically relevant space, which grows and divides. In this thesis, I develop a computationally efficient method for solving reaction-diffusion equations, in a growing and dividing geometry, which approximates growth and division in a rod-shaped bacterium. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-04-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0072720 |
URI | http://hdl.handle.net/2429/42058 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2012-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2012_spring_carlquist_william.pdf [ 2.4MB ]
- Metadata
- JSON: 24-1.0072720.json
- JSON-LD: 24-1.0072720-ld.json
- RDF/XML (Pretty): 24-1.0072720-rdf.xml
- RDF/JSON: 24-1.0072720-rdf.json
- Turtle: 24-1.0072720-turtle.txt
- N-Triples: 24-1.0072720-rdf-ntriples.txt
- Original Record: 24-1.0072720-source.json
- Full Text
- 24-1.0072720-fulltext.txt
- Citation
- 24-1.0072720.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0072720/manifest