A Computationally Efficient Method for Solving Min System Reaction-Diffusion 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 compo- nents. 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 mathemat- ical 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. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 Numerical Solutions to Reaction-Diffusion Equations in Rod- Cell Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1 Previously Published Numerical Approaches for Solving Determinis- tic Min System Reaction-Diffusion Models in vivo . . . . . . . . . . . 5 2.2 Outline of the Base 2 Multi-Time-Step Forward Euler Method and Description within this Chapter . . . . . . . . . . . . . . . . . . . . . 6 2.3 Diffusion in Rod-Cell Geometry . . . . . . . . . . . . . . . . . . . . . 6 2.4 Finite Differences in Rod-Cell Geometry . . . . . . . . . . . . . . . . 8 2.4.1 Deriving Spatial Finite Difference Formulas . . . . . . . . . . 10 2.4.2 Matching and Boundary Conditions . . . . . . . . . . . . . . . 23 2.4.3 The base 2 Multi-Time-Step Forward Euler Method . . . . . . 28 2.4.4 Computational Efficiency of the Base 2 Multi-Time-Step For- ward Euler Method . . . . . . . . . . . . . . . . . . . . . . . . 32 2.4.5 Conservation in the Base 2 Multi-Time-Step Forward Euler Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 iii 2.4.6 Spatial Convergence of the Multi-Time-Step Forward Euler Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.4.7 Temporal Convergence of the Multi-Time-Step Forward Euler Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.4.8 Convergence of the Multi-Time-Step Forward Euler Method to the Infinite domain solution . . . . . . . . . . . . . . . . . . 43 2.4.9 Diffusion on the Rod-Cell Boundary . . . . . . . . . . . . . . . 45 2.5 Reactions in Rod-Cell Geometry . . . . . . . . . . . . . . . . . . . . . 48 3 Numerical Solutions to Reaction-Diffusion Equations on Chang- ing Domains in Rod-Cell Geometry . . . . . . . . . . . . . . . . . 49 3.1 Previously Published Results and Numerical Methods for Solving Reaction-Diffusion Systems in Changing Domains . . . . . . . . . . . 49 3.2 Numerically Solving Reaction-Diffusion Equations in Growing Rod- Cell Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.2.1 Rod-Cell Growth by Adding Cylindrical Sections . . . . . . . 51 3.2.2 Rod-Cell Growth by Stretching the Rod-Cell Body . . . . . . 51 3.2.3 A Comparison of Rod-Cell Growing methods . . . . . . . . . . 52 3.2.4 Hybrid Growth, by Adding Cylindrical Sections and Stretching 52 3.3 Dividing Rod-Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.3.1 The Base 2 Multi-Time-Step Forward Euler Method for Solv- ing Diffusion Equations in Dividing Rod-Cells . . . . . . . . . 54 3.3.2 The Base 2 Multi-Time-Step Forward Euler Method for Solv- ing Diffusion Equations on Dividing Rod-Cell Boundaries . . . 58 3.4 Corrections for Changes in Total Amounts of Conserved Quantities during Rod-Cell Growth and Division . . . . . . . . . . . . . . . . . . 61 4 Examples of Reaction-Diffusion Solutions in Rod-Cell Geometry 62 4.1 The Huang et al. Model . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.2 The Huang et al. Model in a Wild Type Cell that Grows and Divides 64 4.3 The Huang et al. Model in a Continuously Growing Filamentous Mutant Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.4 Future Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 iv List of Tables 2.1 Relative errors for various time-stepping methods based on time step 34 2.2 Base 2 multi-time-step forward Euler method sparsity. . . . . . . . . 35 2.3 Computational time to calculate Matrices required for various methods. 36 2.4 Computational time to calculate one numerical second using various methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 v List of Figures 2.1 Spatial convergence of the base 2 multi-time-step forward Euler method in rod-cell geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.2 Temporal convergence of the base 2 multi-time-step forward Euler method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.3 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 Huang et al. solutions in a wild type cell that grows and divides . . . 65 4.2 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 encour- agement, 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 un- derstandably 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). Ac- cording 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 membrane- bound 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 pole- to-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 (Fischer- Friedrich 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 con- centration 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 trav- eling 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, con- taining 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 Mizu- uchi, 2010). Complex spatial and temporal patterns, arising below the well mixed, constantly flowing solution, demonstrate that spatial and temporal variation in so- lution 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 observa- tions of the Loose et al. experiments (Loose et al., 2008). The model demonstrates traveling and spiral wave solutions, but traveling wave solutions differ 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 buffer 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 behav- iors 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 regu- lar 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 dif- fusion 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 oscilla- tions in 10µm long cells (Huang et al., 2003), behavior consistent with experimen- tally 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. experimen- tal findings (Fischer-Friedrich et al., 2010). Tostevin and Howard proposed a one dimensional stochastic MinD polymerization model, within a dividing domain, pre- dicting that protein diffusion 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 exper- imental 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 pheno- types (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. Nu- merically solving reaction-diffusion models in an in vitro setting is relatively straight forward. However, solving reaction-diffusion equations in a realistic in vivo setting, while maintaining computational efficiency, 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 perti- nent predictions, model behaviors must be determined in three dimensional growing and dividing cells. This thesis outlines the computationally efficient, base 2 multi- time-step forward Euler method, for solving reaction-diffusion equations in growing and dividing rod-cell geometries, and provides examples and applications of numer- ical solutions found using the method. 4 Chapter 2 Numerical Solutions to Reaction-Diffusion Equations in Rod-Cell Geometry Eschericia Coli is a rod-shaped bacterium. To numerically solve reaction-diffusion equations, describing Min system protein dynamics, within an Eschericia Coli in vivo environment, I developed a computationally efficient method, the base 2 multi- time-step forward Euler method, within rod-cell geometry, a geometrical approxi- mation to a rod-shaped bacterium. 2.1 Previously Published Numerical Approaches for Solving Deterministic Min System Reaction- Diffusion Models in vivo Various geometries and numerical methods have been applied to solve deterministic reaction-diffusion 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 differentiation 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 efficient, time-stepping method for solving reaction-diffusion equations, the base 2 multi-time-step forward Euler method. 2.2 Outline of the Base 2 Multi-Time-Step For- ward Euler Method and Description within this Chapter The base 2 multi-time-step forward Euler method for solving reaction-diffusion equa- tions in rod-cell geometry implements multiple time steps, for computationally ef- ficient solution calculations. To solve the diffusion parts of reaction-diffusion equa- tions, detailed in section 2.4, the method utilizes finite difference approximation, described in section 2.4.1, with matching and boundary conditions, defined in sec- tion 2.4.2. Variable time-stepping, within the method, is described in section 2.4.3. Resulting computational efficiency discussion, and comparison of computational ef- ficiency 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 diffusion equations on the rod-cell boundary, which is discussed in section 2.4.9. Finally, complete solutions to reaction-diffusion equations, in rod cell geometry, are calculated by adding changes from reactions to diffusion solutions, in section 2.5. 2.3 Diffusion in Rod-Cell Geometry Reaction-diffusion equations are of type: ∂f(x, t) ∂t = F (f(x, t)) +D∇2f(x, t) , (2.1) 6 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 diffusion constants. To numerically solve systems of type Eq. (2.1), I split each equation within the system into reaction and diffusion parts, then developed a com- putationally efficient method, the base 2 multi-time-step forward Euler method, to numerically solve the diffusion part of each equation within the system: ∂f(x, t) ∂t = D∇2f(x, 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 diffusion 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 diffusion equations, Eq. (2.2), in rod-cell geometry, becomes: ∂fi(r, θ, z, t) ∂t = Di ￿ ∂2fi(r, θ, z, t) ∂r2 + 1 r ∂fi(r, θ, z, t) ∂r + 1 r2 ∂2fi(r, θ, z, t) ∂θ2 + ∂2fi(r, θ, z, t) ∂z2 ￿ , (2.3) within the rod-cell body, ∂gi(r, θ,φ, t) ∂t =Di ￿ ∂2gi(r, θ,φ, t) ∂r2 + 2 r ∂gi(r, θ,φ, t) ∂r + 1 r2 sin2 φ ∂2gi(r, θ,φ, t) ∂θ2 + cosφ r2 sinφ ∂gi(r, θ,φ, t) ∂φ + 1 r2 ∂2gi(r, θ,φ, t) ∂φ2 ￿ , (2.4) within the bottom cap, and ∂hi(r, θ,φ, t) ∂t =Di ￿ ∂2hi(r, θ,φ, t) ∂r2 + 2 r ∂hi(r, θ,φ, t) ∂r + 1 r2 sin2 φ ∂2hi(r, θ,φ, t) ∂θ2 + cosφ r2 sinφ ∂hi(r, θ,φ, t) ∂φ + 1 r2 ∂2hi(r, θ,φ, t) ∂φ2 ￿ , (2.5) 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 repre- sented within both the body and the cap, hence: fi(r, θ, 0) = gi(r, θ , π 2 ), (2.6) 7 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) ∂z = 1 r ∂gi(r, θ , π 2 ) ∂φ , (2.8) ∂2fi(r, θ, 0) ∂z = 1 r ∂gi(r, θ , π 2 ) ∂r + 1 r2 ∂2gi(r, θ , π 2 ) ∂2φ2 , (2.9) Similarly, for the top cap-body interface: ∂fi(r, θ, L) ∂z = 1 r ∂hi(r, θ , π 2 ) ∂φ , (2.10) ∂2fi(r, θ, L) ∂z = 1 r ∂hi(r, θ , π 2 ) ∂r + 1 r2 ∂2hi(r, θ , π 2 ) ∂2φ2 , (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 diffusion 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) ∂r = 0, (2.12) ∂gi(R, θ,φ, t) ∂r = 0, (2.13) ∂hi(R, θ,φ, t) ∂r = 0. (2.14) 2.4 Finite Differences in Rod-Cell Geometry The system of diffusion 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 difference approximation. Finite difference approximation is chosen rather than finite element approximation, as rod-cell geometry is a union of simple geome- tries, which is uniformly discretized fairly easily, and finite difference implementation over various time-steps, for computational efficiency, 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 differ- ences are calculated. Within the rod-cell bottom cap, body, and top cap, there are Nr differing grid points in r, and Nθ differing grid points in θ. Within the rod-cell body, there are Nz differing grid points in z. Within each cap, there are Nφ differing 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 diffusion equations, Eqs. (2.3), (2.4), and (2.5). To avoid singular finite difference values and maintain equal grid point spacing across r = 0, the first grid point in r, of r-index 1, is located at r = ∆r2 , where ∆r is the grid point spacing in r. Hence, ∆r 2 + (Nr − 1)∆r = R. (2.19) ￿ ∆r = R Nr − 12 . (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: ∆θ = 2π Nθ . (2.21) 9 Grid points in z, of z-index l, are ordered from 1, at z = 0, to Nz, at z = L, hence: ∆z = L Nz − 1 . (2.22) 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 diffusion equations, Eqs. (2.3), (2.4), and (2.5), at φ = 0 and φ = π. To avoid singular finite difference 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 φ = ∆φ2 , in the bottom cap, and φ = π − ∆φ2 , in the top cap, where ∆φ is the grid point spacing in φ. Hence, ∆φ 2 + (Nφ − 1)∆φ = π 2 . (2.23) ￿ ∆φ = π 2 Nφ − 12 . (2.24) 2.4.1 Deriving Spatial Finite Difference Formulas Changing variables in diffusion equations from cartesian to cylindrical or spherical geometries results in nonautonomous coefficients, of the form 1r and 1 r2 , in rod-cell body diffusion equation, (2.3), and of the form 2r , 1 r2 sin2 φ , cosφr2 sinφ , and 1 r2 in rod- cell cap diffusion equations, (2.4) and (2.5). These nonautonomous coefficients can increase error in finite difference approximation unless handled carefully. Depending on the value of each nonautonomous coefficient, finite difference 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 difference 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 diffusion equation, (2.3), has a constant coefficient, thus a symmetric three-point finite difference is 10 required for local finite difference approximation error to be of O(∆2): ∂2f(r, θ, z, t) ∂z2 ≈f(r, θ, z −∆z, t) ∆z2 − 2f(r, θ, z, t) ∆z2 + f(r, θ, z +∆z, t) ∆z2 − ∆z 2 12 ∂4f(r, θ, z, t) ∂z4 . (2.25) Nonautonomous coefficients for first order partial derivates of f , g, and h in r, in rod cell diffusion equations, (2.3), (2.4), and (2.5), are of O(r−1). At small r, nonautonomous coefficients of O(r−1) increase error in finite difference approximation by one order of magnitude, O(∆), thus symmetric five point finite differences are required to approximate first order partial derivatives of f , g, and h, in r, with local finite difference approximation errors maximally of O(∆2). Finite differences that approximate second order partial derivatives of f , g, and h, in r, utilize the same grid points as finite differences approximating first order partial derivatives in r, hence: ∂f(r, θ, z, t) ∂r ≈f(r − 2∆r, θ, z, t) 12∆r − 8f(r −∆r, θ, z, t) 12∆r + 8f(r +∆r, θ, z, t) 12∆r − f(r + 2∆r, θ, z, t) 12∆r + ∆r4 30 ∂5f(r, θ, z, t) ∂r5 , (2.26) ∂2f(r, θ, z, t) ∂r2 ≈−f(r − 2∆r, θ, z, t) 12∆r2 + 16f(r −∆r, θ, z, t) 12∆r2 − 30f(r, θ, z, t) 12∆r2 + 16f(r +∆r, θ, z, t) 12∆r2 − f(r + 2∆r, θ, z, t) 12∆r2 + ∆r4 90 ∂5f(r, θ, z, t) ∂r6 , (2.27) ∂g(r, θ,φ, t) ∂r ≈g(r − 2∆r, θ,φ, t) 12∆r − 8g(r −∆r, θ,φ, t) 12∆r + 8g(r +∆r, θ,φ, t) 12∆r − g(r + 2∆r, θ,φ, t) 12∆r + ∆r4 30 ∂5g(r, θ,φ, t) ∂r5 , (2.28) ∂2g(r, θ,φ, t) ∂r2 ≈−g(r − 2∆r, θ,φ, t) 12∆r2 + 16g(r −∆r, θ,φ, t) 12∆r2 − 30g(r, θ,φ, t) 12∆r2 + 16g(r +∆r, θ,φ, t) 12∆r2 − g(r + 2∆r, θ,φ, t) 12∆r2 + ∆r4 90 ∂5g(r, θ,φ, t) ∂r6 , (2.29) 11 ∂h(r, θ,φ, t) ∂r ≈h(r − 2∆r, θ,φ, t) 12∆r − 8h(r −∆r, θ,φ, t) 12∆r + 8h(r +∆r, θ,φ, t) 12∆r − h(r + 2∆r, θ,φ, t) 12∆r + ∆r4 30 ∂5h(r, θ,φ, t) ∂r5 , (2.30) ∂2h(r, θ,φ, t) ∂r2 ≈−h(r − 2∆r, θ,φ, t) 12∆r2 + 16h(r −∆r, θ,φ, t) 12∆r2 − 30h(r, θ,φ, t) 12∆r2 + 16h(r +∆r, θ,φ, t) 12∆r2 − h(r + 2∆r, θ,φ, t) 12∆r2 + ∆r4 90 ∂5h(r, θ,φ, t) ∂r6 . (2.31) Within the rod-cell body, grid points in r are located directly across the r origin, differing by θ = π, thus may act as grid points in the negative r direction, for finite difference 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 difference formulas that require no grid points in the negative r direction. For r index 1: ∂g(r, θ,φ, t) ∂r ≈−25g(r, θ,φ, t) 12∆r + 48g(r +∆r, θ,φ, t) 12∆r − 36g(r + 2∆r, θ,φ, t) 12∆r + 16g(r + 3∆r, θ,φ, t) 12∆r − 3g(r + 4∆r, θ,φ, t) 12∆r + ∆r4 5 ∂5g(r, θ,φ, t) ∂r5 , (2.32) ∂2g(r, θ,φ, t) ∂r2 ≈35g(r, θ,φ, t) 12∆r2 − 104g(r +∆r, θ,φ, t) 12∆r2 + 114g(r + 2∆r, θ,φ, t) 12∆r2 − 56g(r + 3∆r, θ,φ, t) 12∆r2 + 11g(r + 4∆r, θ,φ, t) 12∆r2 + 5∆r4 6 ∂5g(r, θ,φ, t) ∂r5 , (2.33) ∂h(r, θ,φ, t) ∂r ≈−25h(r, θ,φ, t) 12∆r + 48h(r +∆r, θ,φ, t) 12∆r − 36h(r + 2∆r, θ,φ, t) 12∆r + 16h(r + 3∆r, θ,φ, t) 12∆r − 3h(r + 4∆r, θ,φ, t) 12∆r + ∆r4 5 ∂5h(r, θ,φ, t) ∂r5 , (2.34) ∂2h(r, θ,φ, t) ∂r2 ≈35h(r, θ,φ, t) 12∆r2 − 104h(r +∆r, θ,φ, t) 12∆r2 + 114h(r + 2∆r, θ,φ, t) 12∆r2 − 56h(r + 3∆r, θ,φ, t) 12∆r2 + 11h(r + 4∆r, θ,φ, t) 12∆r2 − 5∆r 4 6 ∂5h(r, θ,φ, t) ∂r5 . (2.35) 12 For r index 2: ∂g(r, θ,φ, t) ∂r ≈−3g(r −∆r, θ,φ, t) 12∆r − 10g(r, θ,φ, t) 12∆r + 18g(r +∆r, θ,φ, t) 12∆r − 6g(r + 2∆r, θ,φ, t) 12∆r + g(r + 3∆r, θ,φ, t) 12∆r + ∆r4 20 ∂5g(r, θ,φ, t) ∂r5 , (2.36) ∂2g(r, θ,φ, t) ∂r2 ≈11g(r −∆r, θ,φ, t) 12∆r2 − 20g(r, θ,φ, t) 12∆r2 + 6g(r +∆r, θ,φ, t) 12∆r2 + 4g(r + 2∆r, θ,φ, t) 12∆r2 − g(r + 3∆r, θ,φ, t) 12∆r2 − ∆r 3 12 ∂5g(r, θ,φ, t) ∂r5 , (2.37) ∂h(r, θ,φ, t) ∂r ≈−3h(r −∆r, θ,φ, t) 12∆r − 10h(r, θ,φ, t) 12∆r + 18h(r +∆r, θ,φ, t) 12∆r − 6h(r + 2∆r, θ,φ, t) 12∆r + h(r + 3∆r, θ,φ, t) 12∆r + ∆r4 20 ∂5h(r, θ,φ, t) ∂r5 , (2.38) ∂2h(r, θ,φ, t) ∂r2 ≈11h(r −∆r, θ,φ, t) 12∆r2 − 20h(r, θ,φ, t) 12∆r2 + 6h(r +∆r, θ,φ, t) 12∆r2 + 4h(r + 2∆r, θ,φ, t) 12∆r2 − h(r + 3∆r, θ,φ, t) 12∆r2 − ∆r 3 12 ∂5h(r, θ,φ, t) ∂r5 . (2.39) The nonautonomous coefficient for the second order partial derivative of f in θ, in rod cell diffusion equation, (2.3), is of O(r−2), which, for small r, increase local error from finite difference approximation by two orders of magnitude, O(∆2), thus a symmetric five-point finite difference is required to approximate the second order partial derivative of f in θ with local finite difference approximation error of O(∆2): ∂2f(r, θ, z, t) ∂θ2 ≈−f(r, θ − 2∆θ, z, t) 12∆θ2 + 16f(r, θ −∆θ, z, t) 12∆θ2 − 30f(r, θ, z, t) 12∆θ2 + 16f(r, θ +∆θ, z, t) 12∆θ2 − f(r, θ + 2∆θ, z, t) 12∆θ2 + ∆θ4 90 ∂5f(r, θ, z, t) ∂θ6 . (2.40) Nonautonomous coefficients for second order partial derivative of g and h in θ, in rod cell diffusion equations (2.4) and (2.5), r−2 sin−2(φ), increases local 13 error from finite difference approximation, at small r and small φ, by four orders of magnitude, O(∆4), thus symmetric seven-point finite differences are required to approximate second order partial derivatives of g and h in θ with local finite difference approximation errors maximally of O(∆2): ∂2g(r, θ,φ, t) ∂θ2 ≈2g(r, θ − 3∆θ,φ, t) 180∆θ2 − 27g(r, θ − 2∆θ,φ, t) 180∆θ2 + 270g(r, θ −∆θ,φ, t) 180∆θ2 − 490g(r, θ,φ, t) 180∆θ2 + 270g(r, θ +∆θ,φ, t) 180∆θ2 − 27g(r, θ + 2∆θ,φ, t) 180∆θ2 + 2g(r, θ + 3∆θ,φ, t) 180∆θ2 − ∆θ 6 560 ∂8g(r, θ,φ, t) ∂θ8 . (2.41) ∂2h(r, θ,φ, t) ∂θ2 ≈2h(r, θ − 3∆θ,φ, t) 180∆θ2 − 27h(r, θ − 2∆θ,φ, t) 180∆θ2 + 270h(r, θ −∆θ,φ, t) 180∆θ2 − 490h(r, θ,φ, t) 180∆θ2 + 270h(r, θ +∆θ,φ, t) 180∆θ2 − 27h(r, θ + 2∆θ,φ, t) 180∆θ2 + 2h(r, θ + 3∆θ,φ, t) 180∆θ2 − ∆θ 6 560 ∂8h(r, θ,φ, t) ∂θ8 . (2.42) Partial derivatives of g and h in φ, in rod cell diffusion equations (2.4) and (2.5), have nonautonomous coefficients cos(φ)r2 sin(φ) and r −2. At small r and small φ, the nonautonomous coefficient cos(φ)r2 sin(φ) increases local errors from finite difference approximation by three orders of magnitude, thus symmetric seven-point finite dif- ferences of g and h in φ are required for local finite difference 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 inter- face, thus symmetric seven-point finite differences that approximate first and second 14 order partial derivatives of g and h in φ are implemented: ∂g(r, θ,φ, t) ∂φ ≈−g(r, θ,φ− 3∆φ, t) 60∆φ2 + 9g(r, θ,φ− 2∆φ, t) 60∆φ − 45g(r, θ,φ−∆φ, t) 60∆φ + 45g(r, θ,φ+∆φ, t) 60∆φ − 9g(r, θ,φ+ 2∆φ, t) 60∆φ + g(r, θ,φ+ 3∆φ, t) 60∆φ − ∆φ 6 140 ∂7g(r, θ,φ, t) ∂φ7 , (2.43) ∂2g(r, θ,φ, t) ∂φ2 ≈2g(r, θ,φ− 3∆φ, t) 180∆φ2 − 27g(r, θ,φ− 2∆φ, t) 180∆φ2 + 270g(r, θ,φ−∆φ, t) 180∆φ2 − 490g(r, θ,φ, t) 180∆φ2 + 270g(r, θ,φ+∆φ, t) 180∆φ2 − 27g(r, θ,φ+ 2∆φ, t) 180∆φ2 + 2g(r, θ,φ+ 3∆φ, t) 180∆φ2 − ∆φ 6 560 ∂8g(r, θ,φ, t) ∂φ8 , (2.44) ∂h(r, θ,φ, t) ∂φ ≈−h(r, θ,φ− 3∆φ, t) 60∆φ2 + 9h(r, θ,φ− 2∆φ, t) 60∆φ − 45h(r, θ,φ−∆φ, t) 60∆φ + 45h(r, θ,φ+∆φ, t) 60∆φ − 9h(r, θ,φ+ 2∆φ, t) 60∆φ + h(r, θ,φ+ 3∆φ, t) 60∆φ − ∆φ 6 140 ∂7h(r, θ,φ, t) ∂φ7 , (2.45) ∂2h(r, θ,φ, t) ∂φ2 ≈2h(r, θ,φ− 3∆φ, t) 180∆φ2 − 27h(r, θ,φ− 2∆φ, t) 180∆φ2 + 270h(r, θ,φ−∆φ, t) 180∆φ2 − 490h(r, θ,φ, t) 180∆φ2 + 270h(r, θ,φ+∆φ, t) 180∆φ2 − 27h(r, θ,φ+ 2∆φ, t) 180∆φ2 + 2h(r, θ,φ+ 3∆φ, t) 180∆φ2 − ∆φ 6 560 ∂8h(r, θ,φ, t) ∂φ8 , (2.46) Within each cap, at grid points of φ-index, Nφ− 2, the nonautonomous coef- ficient, cos(φ)r2 sin(φ , is of O(r −2), thus symmetric five-point finite differences that approx- imate first and second order partial derivatives of g and h in φ, in rod cell diffusion equations (2.4) and (2.5), are required for local finite difference approximation errors 15 to be maximally of O(∆2): ∂g(r, θ,φ, t) ∂φ ≈g(r, θ,φ− 2∆φ, t) 12∆φ − 8g(r, θ,φ−∆φ, t) 12∆φ + 8g(r, θ,φ+∆φ, t) 12∆φ − g(r, θ,φ+ 2∆φ, t) 12∆φ + ∆φ4 30 ∂5g(r, θ,φ, t) ∂φ5 , (2.47) ∂2g(r, θ,φ, t) ∂φ2 ≈−g(r, θ,φ− 2∆φ, t) 12∆φ2 + 16g(r, θ,φ−∆φ, t) 12∆φ2 − 30g(r, θ,φ, t) 12∆φ2 + 16g(r, θ,φ+∆φ, t) 12∆φ2 − g(r, θ,φ+ 2∆φ, t) 12∆φ2 + ∆φ4 90 ∂6g(r, θ,φ, t) ∂φ6 , (2.48) ∂h(r, θ,φ, t) ∂φ ≈h(r, θ,φ− 2∆φ, t) 12∆φ − 8h(r, θ,φ−∆φ, t) 12∆φ + 8h(r, θ,φ+∆φ, t) 12∆φ − h(r, θ,φ+ 2∆φ, t) 12∆φ + ∆φ4 30 ∂5h(r, θ,φ, t) ∂φ5 , (2.49) ∂2h(r, θ,φ, t) ∂φ2 ≈−h(r, θ,φ− 2∆φ, t) 12∆φ2 + 16h(r, θ,φ−∆φ, t) 12∆φ2 − 30h(r, θ,φ, t) 12∆φ2 + 16h(r, θ,φ+∆φ, t) 12∆φ2 − h(r, θ,φ+ 2∆φ, t) 12∆φ2 + ∆φ4 90 ∂6h(r, θ,φ, t) ∂φ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 difference function value across the interface, for finite differences in z and φ (derived in Section 2.4.2). Thus, finite differences in z and φ are calculated with at most one finite difference function value across the interface. Within each cap, at grid points of φ-index, Nφ − 1, the nonlinear partial derivative multiplier, cos(φ)r2 sin(φ) , is of O(r−2). Implementing five point finite differences of g and h in φ, Eqs 2.47, 2.48, 2.49, and 2.50, forces local finite difference approximation errors to be maximally of O(∆2), and require only one finite difference function value across each interface. Within each cap, at grid points within the cap-body interface, of φ-index, Nφ, the nonautonomous coefficients for the second order partial derivatives of g and h in φ, in rod cell diffusion equations (2.4) and (2.5), are of O(r−2). Finite differences may only implement one function value across each interface, calculated from interface matching conditions. Thus, non-symmetric six-point finite differences are required to approximate first and second order partial derivatives of g and h in 16 φ, with local finite difference approximation errors maximally of O(∆2): ∂g(r, θ,φ, t) ∂φ ≈3g(r, θ,φ− 4∆φ, t) 60∆φ2 − 20g(r, θ,φ− 3∆φ, t) 60∆φ2 + 60g(r, θ,φ− 2∆φ, t) 60∆φ − 120g(r, θ,φ−∆φ, t) 60∆φ + 65g(r, θ,φ, t) 60∆φ + 12g(r, θ,φ+∆φ, t) 60∆φ − ∆φ 5 30 ∂6g(r, θ,φ, t) ∂φ6 , (2.51) ∂2g(r, θ,φ, t) ∂φ2 ≈g(r, θ,φ− 4∆φ, t) 12∆φ2 − 6g(r, θ,φ− 3∆φ, t) 12∆φ2 + 14g(r, θ,φ− 2∆φ, t) 12∆φ2 − 4g(r, θ,φ−∆φ, t) 12∆φ2 − 15g(r, θ,φ, t) 12∆φ2 + 10g(r, θ,φ+∆φ, t) 12∆φ2 − 13∆φ 4 180 ∂6g(r, θ,φ, t) ∂φ6 , (2.52) ∂h(r, θ,φ, t) ∂φ ≈3h(r, θ,φ− 4∆φ, t) 60∆φ2 − 20h(r, θ,φ− 3∆φ, t) 60∆φ2 + 60h(r, θ,φ− 2∆φ, t) 60∆φ − 120h(r, θ,φ−∆φ, t) 60∆φ + 65h(r, θ,φ, t) 60∆φ + 12h(r, θ,φ+∆φ, t) 60∆φ − ∆φ 5 30 ∂6h(r, θ,φ, t) ∂φ6 , (2.53) ∂2h(r, θ,φ, t) ∂φ2 ≈h(r, θ,φ− 4∆φ, t) 12∆φ2 − 6h(r, θ,φ− 3∆φ, t) 12∆φ2 + 14h(r, θ,φ− 2∆φ, t) 12∆φ2 − 4h(r, θ,φ−∆φ, t) 12∆φ2 − 15h(r, θ,φ, t) 12∆φ2 + 10h(r, θ,φ+∆φ, t) 12∆φ2 − 13∆φ 4 180 ∂6h(r, θ,φ, t) ∂φ6 , (2.54) Finite Difference 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: ∆2rfi(rj, θk, zl, tn) = −fi(rj−2, θk, zl, tn) 12∆r2 + 16fi(rj−1, θk, zl, tn) 12∆r2 − 30fi(rj, θk, zl, tn) 12∆r2 + 16fi(rj+1, θk, zl, tn) 12∆r2 − fi(rj+2, θk, zl, tn) 12∆r2 , (2.55) 17 for j ∈ {1, . . . , Nr}. ∆1rfi(rj, θk, zl, tn) = fi(rj−2, θk, zl, tn) 12∆r − 8fi(rj−1, θk, zl, tn) 12∆r + 8fi(rj+1, θk, zl, tn) 12∆r (2.56) − fi(rj+2, θk, zl, tn) 12∆r , (2.57) for j ∈ {1, . . . , Nr}. ∆2θfi(rj, θk, zl, tn) = fi(rj, θk−2, zl, tn) 12∆θ2 + 16fi(rj, θk−1, zl, tn) 12∆θ2 − 30fi(rj, θk, zl, tn) 12∆θ2 + 16fi(rj, θk+1, zl, tn) 12∆θ2 − fi(rj, θk+2, zl, tn) 12∆θ2 , (2.58) for k ∈ {1, . . . , Nθ}. ∆2zfi(rj, θk, zl, tn) = fi(rj, θk, zl−1, tn) ∆z2 − 2fi(rj, θk, zl, tn) ∆z2 + fi(rj, θk, zl+1, tn) ∆z2 , (2.59) for l ∈ {1, . . . , Nz}. ∆2rgi(rj, θk,φm, tn) = 35gi(rj, θk,φm, tn) 12∆r2 − 104gi(rj+1, θk,φm, tn) 12∆r2 + 114gi(rj+2, θk,φm, tn) 12∆r2 − 56gi(rj+3, θk,φm, tn) 12∆r2 + 11gi(rj+4, θk,φm, tn) 12∆r2 , (2.60) for j ∈ {1}. ∆2rgi(rj, θk,φm, tn) = 11gi(rj−1, θk,φm, tn) 12∆r2 − 20gi(rj, θk,φm, tn) 12∆r2 + 6gi(rj+1, θk,φm, tn) 12∆r2 + 4gi(rj+2, θk,φm, tn) 12∆r2 − gi(rj+3, θk,φm, tn) 12∆r2 , (2.61) for j ∈ {2}. ∆2rgi(rj, θk,φm, tn) = −gi(rj−2, θk,φm, tn) 12∆r2 + 16gi(rj−1, θk,φm, tn) 12∆r2 − 30gi(rj, θk,φm, tn) 12∆r2 + 16gi(rj+1, θk,φm, tn) 12∆r2 − gi(rj+2, θk,φm, tn) 12∆r2 , (2.62) for j ∈ {3, . . . , Nr}. ∆1rgi(rj, θk,φm, tn) = −25gi(rj, θk,φm, tn) 12∆r + 48gi(rj+1, θk,φm, tn) 12∆r − 36gi(rj+2, θk,φm, tn) 12∆r + 16gi(rj+3, θk,φm, tn) 12∆r − 3gi(rj+4, θk,φm, tn) 12∆r , (2.63) 18 for j ∈ {1}. ∆1rgi(rj, θk,φm, tn) = −3gi(rj−1, θk,φm, tn) 12∆r − 10gi(rj, θk,φm, tn) 12∆r + 18gi(rj+1, θk,φm, tn) 12∆r − 6gi(rj+2, θk,φm, tn) 12∆r + gi(rj+3, θk,φm, tn) 12∆r , (2.64) for j ∈ {2}. ∆1rgi(rj, θk,φm, tn) = gi(rj−2, θk,φm, tn) 12∆r − 8gi(rj−1, θk,φm, tn) 12∆r + 8gi(rj+1, θk,φm, tn) 12∆r (2.65) − gi(rj+2, θk,φm, tn) 12∆r , (2.66) for j ∈ {3, . . . , Nr}. ∆2θgi(rj, θk,φm, tn) = 2gi(rj, θk−3,φm, tn) 180∆θ2 − 27gi(rj, θk−2,φm, tn) 180∆θ2 + 270gi(rj, θk−1,φm, tn) 180∆θ2 − 490gi(rj, θk,φm, tn) 180∆θ2 + 270gi(rj, θk+1,φm, tn) 180∆θ2 − 27gi(rj, θk+2,φm, tn) 180∆θ2 + 2gi(rj, θk+3,φm, tn) 180∆θ2 , (2.67) for k ∈ {1, . . . , Nθ}. ∆1φgi(rj, θk,φm, tn) = −gi(rj, θk,φm−3, tn) 60∆φ + 9gi(rj, θk,φm−2, tn) 60∆φ − 45gi(rj, θk,φm−1, tn) 60∆φ − 45gi(rj, θk,φm+1, tn) 60∆φ − 9gi(rj, θk,φm+2, tn) 60∆φ + gi(rj, θk,φm+3, tn) 60∆φ , (2.68) for m ∈ {1, . . . , Nφ − 3}. ∆1φgi(rj, θk,φm, tn) = gi(rj, θk,φm−2, tn) 12∆φ − 8gi(rj, θk,φm−1, tn) 12∆φ + 8gi(rj, θk,φm+1, tn) 12∆φ − gi(rj, θk,φm+2, tn) 12∆φ , (2.69) for m ∈ {Nφ − 2, Nφ − 1}. ∆1φgi(rj, θk,φm, tn) = 3gi(rj, θk,φm−4, tn) 60∆φ − 20gi(rj, θk,φm−3, tn) 60∆φ + 60gi(rj, θk,φm−2, tn) 60∆φ − 120gi(rj, θk,φm−1, tn) 60∆φ + 65gi(rj, θk,φm, tn) 60∆φ + 12gi(rj, θk,φm+1, tn) 60∆φ , (2.70) 19 for m ∈ {Nφ}. ∆2φgi(rj, θk,φm, tn) = 2gi(rj, θk,φm−3, tn) 180∆φ2 − 27gi(rj, θk,φm−2, tn) 180∆φ2 + 270gi(rj, θk,φm−1, tn) 180∆φ2 − 490gi(rj, θk,φm, tn) 180∆φ2 + 270gi(rj, θk,φm+1, tn) 180∆φ2 − 27gi(rj, θk,φm+2, tn) 180∆φ2 + 2gi(rj, θk,φm+3, tn) 180∆φ2 , (2.71) for m ∈ {1, . . . , Nφ − 3}. ∆2φgi(rj, θk,φm, tn) = −gi(rj, θk,φm−2, tn) 12∆φ2 + 16gi(rj, θk,φm−1, tn) 12∆φ2 − 30gi(rj, θk,φm, tn) 12∆φ2 + 16gi(rj, θk,φm+1, tn) 12∆φ2 − gi(rj, θk,φm+2, tn) 12∆φ2 , (2.72) for m ∈ {Nφ − 2, Nφ − 1}. ∆2φgi(rj, θk,φm, tn) = gi(rj, θk,φm−4, tn) 12∆φ2 − 6gi(rj, θk,φm−3, tn) 12∆φ2 + 14gi(rj, θk,φm−2, tn) 12∆φ2 − 4gi(rj, θk,φm−1, tn) 12∆φ2 − 15gi(rj, θk,φm, tn) 12∆φ2 + 10gi(rj, θk,φm+1, tn) 12∆φ2 , (2.73) for m ∈ {Nφ}. ∆2rhi(rj, θk,φm, tn) = 35hi(rj, θk,φm, tn) 12∆r2 − 104hi(rj+1, θk,φm, tn) 12∆r2 + 114hi(rj+2, θk,φm, tn) 12∆r2 − 56hi(rj+3, θk,φm, tn) 12∆r2 + 11hi(rj+4, θk,φm, tn) 12∆r2 , (2.74) for j ∈ {1}. ∆2rhi(rj, θk,φm, tn) = 11hi(rj−1, θk,φm, tn) 12∆r2 − 20hi(rj, θk,φm, tn) 12∆r2 + 6hi(rj+1, θk,φm, tn) 12∆r2 + 4hi(rj+2, θk,φm, tn) 12∆r2 − hi(rj+3, θk,φm, tn) 12∆r2 , (2.75) for j ∈ {2}. ∆2rhi(rj, θk,φm, tn) = −hi(rj−2, θk,φm, tn) 12∆r2 + 16hi(rj−1, θk,φm, tn) 12∆r2 − 30hi(rj, θk,φm, tn) 12∆r2 + 16hi(rj+1, θk,φm, tn) 12∆r2 − hi(rj+2, θk,φm, tn) 12∆r2 , (2.76) 20 for j ∈ {3, . . . , Nr}. ∆1rhi(rj, θk,φm, tn) = −25hi(rj, θk,φm, tn) 12∆r + 48hi(rj+1, θk,φm, tn) 12∆r − 36hi(rj+2, θk,φm, tn) 12∆r + 16hi(rj+3, θk,φm, tn) 12∆r − 3hi(rj+4, θk,φm, tn) 12∆r , (2.77) for j ∈ {1}. ∆1rhi(rj, θk,φm, tn) = −3hi(rj−1, θk,φm, tn) 12∆r − 10hi(rj, θk,φm, tn) 12∆r + 18hi(rj+1, θk,φm, tn) 12∆r − 6hi(rj+2, θk,φm, tn) 12∆r + hi(rj+3, θk,φm, tn) 12∆r , (2.78) for j ∈ {2}. ∆1rhi(rj, θk,φm, tn) = hi(rj−2, θk,φm, tn) 12∆r − 8hi(rj−1, θk,φm, tn) 12∆r + 8hi(rj+1, θk,φm, tn) 12∆r (2.79) − hi(rj+2, θk,φm, tn) 12∆r , (2.80) for j ∈ {3, . . . , Nr}. ∆2θhi(rj, θk,φm, tn) = 2hi(rj, θk−3,φm, tn) 180∆θ2 − 27hi(rj, θk−2,φm, tn) 180∆θ2 + 270hi(rj, θk−1,φm, tn) 180∆θ2 − 490hi(rj, θk,φm, tn) 180∆θ2 + 270hi(rj, θk+1,φm, tn) 180∆θ2 − 27hi(rj, θk+2,φm, tn) 180∆θ2 + 2hi(rj, θk+3,φm, tn) 180∆θ2 , (2.81) for k ∈ {1, . . . , Nθ}. ∆1φhi(rj, θk,φm, tn) = −hi(rj, θk,φm−3, tn) 60∆φ + 9hi(rj, θk,φm−2, tn) 60∆φ − 45hi(rj, θk,φm−1, tn) 60∆φ − 45hi(rj, θk,φm+1, tn) 60∆φ − 9hi(rj, θk,φm+2, tn) 60∆φ + hi(rj, θk,φm+3, tn) 60∆φ , (2.82) for m ∈ {1, . . . , Nφ − 3}. ∆1φhi(rj, θk,φm, tn) = hi(rj, θk,φm−2, tn) 12∆φ − 8hi(rj, θk,φm−1, tn) 12∆φ + 8hi(rj, θk,φm+1, tn) 12∆φ − hi(rj, θk,φm+2, tn) 12∆φ , (2.83) 21 for m ∈ {Nφ − 2, Nφ − 1}. ∆1φhi(rj, θk,φm, tn) = 3hi(rj, θk,φm−4, tn) 60∆φ − 20hi(rj, θk,φm−3, tn) 60∆φ + 60hi(rj, θk,φm−2, tn) 60∆φ − 120hi(rj, θk,φm−1, tn) 60∆φ + 65hi(rj, θk,φm, tn) 60∆φ + 12hi(rj, θk,φm+1, tn) 60∆φ , (2.84) for m ∈ {Nφ}. ∆2φhi(rj, θk,φm, tn) = 2hi(rj, θk,φm−3, tn) 180∆φ2 − 27hi(rj, θk,φm−2, tn) 180∆φ2 + 270hi(rj, θk,φm−1, tn) 180∆φ2 − 490hi(rj, θk,φm, tn) 180∆φ2 + 270hi(rj, θk,φm+1, tn) 180∆φ2 − 27hi(rj, θk,φm+2, tn) 180∆φ2 + 2hi(rj, θk,φm+3, tn) 180∆φ2 , (2.85) for m ∈ {1, . . . , Nφ − 3}. ∆2φhi(rj, θk,φm, tn) = −hi(rj, θk,φm−2, tn) 12∆φ2 + 16hi(rj, θk,φm−1, tn) 12∆φ2 − 30hi(rj, θk,φm, tn) 12∆φ2 + 16hi(rj, θk,φm+1, tn) 12∆φ2 − hi(rj, θk,φm+2, tn) 12∆φ2 , (2.86) for m ∈ {Nφ − 2, Nφ − 1}. ∆2φhi(rj, θk,φm, tn) = hi(rj, θk,φm−4, tn) 12∆φ2 − 6hi(rj, θk,φm−3, tn) 12∆φ2 + 14hi(rj, θk,φm−2, tn) 12∆φ2 − 4hi(rj, θk,φm−1, tn) 12∆φ2 − 15hi(rj, θk,φm, tn) 12∆φ2 + 10hi(rj, θk,φm+1, tn) 12∆φ2 , (2.87) for m ∈ {Nφ}. Spatial Finite Differences Grid points at the interface of the rod-cell body and either cap are represented in both cylindrical and spherical coordinates, thus diffusion occurs within both cylindrical and spherical geometry, and is invariant from matching conditions. To 22 ensure that diffusivity is not double the intended diffusivity, the total diffusion, at the interface of the rod-cell body and either cap, is half cylindrical diffusion and half spherical diffusion. Spatial finite difference formulas for the system of diffusion 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 ≈ Di ￿￿ 1− δl1 + δlNz 2 ￿ · ￿ ∆2rfi(rj, θk, zl, tn) + 1 rj ∆1rfi(rj, θk, zl, tn) + 1 r2j ∆2θfi(rj, θk, zl, tn) +∆ 2 zfi(rj, θk, zl, tn) ￿￿ (2.88) ∂gi(rj, θk,φm, tn) ∂t ≈ Di ￿￿ 1− δmNφ 2 ￿ · ￿ ∆2rgi(rj, θk,φm, tn) + 2 rj ∆1rgi(rj, θk,φm, tn) + 1 r2j sin 2 φm ∆2θgi(rj, θk,φm, tn) + cosφm r2j sinφm ∆1φgi(rj, θk,φm, tn) + 1 r2j ∆2φgi(rj, θk,φm, tn) ￿￿ (2.89) ∂hi(rj, θk,φm, tn) ∂t ≈ Di ￿￿ 1− δmNφ 2 ￿ · ￿ ∆2rhi(rj, θk,φm, tn) + 2 rj ∆1rhi(rj, θk,φm, tn) + 1 r2j sin 2 φm ∆2θhi(rj, θk,φm, tn) + cosφm r2j sinφm ∆1φhi(rj, θk,φm, tn) + 1 r2j ∆2φhi(rj, θk,φm, tn) ￿￿ (2.90) where δxy is the Kronecker delta, with spatial variable indices x and y. 2.4.2 Matching and Boundary Conditions Finite difference boundary conditions are chosen to match those described in sec- tion 2.3, for diffusion 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 difference approximations of bottom cap-body matching conditions, 2.8 and 2.9, matching first and second order partial derivative in z, implementing the ap- propriate number of grid points for approximations to be maximally of O(∆2): 0 = fi(rj, θk, z0, tn) 2∆z − fi(rj, θk, z2, tn) 2∆z + 3gi(rj, θk,φNφ−4, tn) 60rj∆φ − 20gi(rj, θk,φNφ−3, tn) 60rj∆φ + 60gi(rj, θk,φNφ−2, tn) 60rj∆φ − 120gi(rj, θk,φNφ−1, tn) 60rj∆φ + 65gi(rj, θk,φNφ , tn) 60rj∆φ + 12gi(rj, θk,φNφ+1, tn) 60rj∆φ , (2.93) and 0 = −fi(rj, θk, z0, tn) ∆z2 + 2fi(rj, θk, z1, tn) ∆z2 − fi(rj, θk, z2, tn) ∆z2 + gi(rj−2, θk,φNφ , tn) 12rj∆r − 8gi(rj−1, θk,φNφ , tn) 12rj∆r + 8gi(rj+1, θk,φNφ , tn) 12rj∆r − gi(rj+2, θk,φNφ , tn) 12rj∆r + gi(rj, θk,φNφ−4, tn) 12r2j∆φ 2 − 6gi(rj, θk,φNφ−3, tn) 12r2j∆φ 2 + 14gi(rj, θk,φNφ−2, tn) 12r2j∆φ 2 − 4gi(rj, θk,φNφ−1, tn) 12r2j∆φ 2 − 15gi(rj, θk,φNφ , tn) 12r2j∆φ 2 + 10gi(rj, θk,φNφ+1, tn) 12r2j∆φ 2 , (2.94) 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 difference formulas. fi(rj,θk, z0, tn) = 24rj∆φ 25∆z + 12rj∆φ fi(rj, θk, z1, tn) + ￿ 50∆z 25∆z + 12rj∆φ + 1 ￿ fi(rj, θk, z2, tn) − 3∆z 2 2rj(25∆z + 12rj∆φ) gi(rj, θk,φNφ−4, tn) + 32∆z2 3rj(25∆z + 12rj∆φ) gi(rj, θk,φNφ−3, tn) − 36∆z 2 rj(25∆z + 12rj∆φ) gi(rj, θk,φNφ−2, tn) + 96∆z2 rj(25∆z + 12rj∆φ) gi(rj, θk,φNφ−1, tn) − 415∆z 2 6rj(25∆z + 12rj∆φ) gi(rj, θk,φNφ , tn) + ∆φ∆z2 ∆r(25∆z + 12rj∆φ) gi(rj−2, θk,φNφ , tn) − 8∆φ∆z 2 ∆r(25∆z + 12rj∆φ) gi(rj−1, θk,φNφ , tn) + 8∆φ∆z2 ∆r(25∆z + 12rj∆φ) gi(rj+1, θk,φNφ , tn) − ∆φ∆z 2 ∆r(25∆z + 12rj∆φ) gi(rj+2, θk,φNφ , tn) (2.95) gi(rj, θk,φNφ+1, tn) = −60r2j∆φ2 ∆z(25∆z + 12rj∆φ) fi(rj, θk, z1, tn) + 60r2j∆φ 2 ∆z(25∆z + 12rj∆φ) fi(rj, θk, z2, tn) + ￿ 15∆z 4(25∆z + 12rj∆φ) − 1 4 ￿ gi(rj, θk,φNφ−4, tn) − ￿ 80∆z 3(25∆z + 12rj∆φ) + 5 3 ￿ gi(rj, θk,φNφ−3, tn) + ￿ 90∆z 25∆z + 12rj∆φ − 5 ￿ gi(rj, θk,φNφ−2, tn) − ￿ 240∆z 25∆z + 12rj∆φ + 10 ￿ gi(rj, θk,φNφ−1, tn) + ￿ 2075∆z 12(25∆z + 12rj∆φ) − 65 12 ￿ gi(rj, θk,φNφ , tn) − 5rj∆φ 2∆z 2∆r(25∆z + 12rj∆φ) gi(rj−2, θk,φNφ , tn) + 20rj∆φ2∆z ∆r(25∆z + 12rj∆φ) gi(rj−1, θk,φNφ , tn) − 20rj∆φ 2∆z ∆r(25∆z + 12rj∆φ) gi(rj+1, θk,φNφ , tn) + 5rj∆φ2∆z 2∆r(25∆z + 12rj∆φ) gi(rj+2, θk,φNφ , tn) (2.96) Similarly, by combining finite difference 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: fi(rj,θk, zNz+1, tn) = 24rj∆φ 25∆z + 12rj∆φ fi(rj, θk, zNz , tn) + ￿ 50∆z 25∆z + 12rj∆φ + 1 ￿ fi(rj, θk, zNz−1, tn)− 3∆z2 2rj(25∆z + 12rj∆φ) hi(rj, θk,φNφ−4, tn) + 32∆z2 3rj(25∆z + 12rj∆φ) hi(rj, θk,φNφ−3, tn)− 36∆z2 rj(25∆z + 12rj∆φ) hi(rj, θk,φNφ−2, tn) + 96∆z2 rj(25∆z + 12rj∆φ) hi(rj, θk,φNφ−1, tn)− 415∆z2 6rj(25∆z + 12rj∆φ) hi(rj, θk,φNφ , tn) + ∆φ∆z2 ∆r(25∆z + 12rj∆φ) hi(rj−2, θk,φNφ , tn)− 8∆φ∆z2 ∆r(25∆z + 12rj∆φ) hi(rj−1, θk,φNφ , tn) + 8∆φ∆z2 ∆r(25∆z + 12rj∆φ) hi(rj+1, θk,φNφ , tn)− ∆φ∆z2 ∆r(25∆z + 12rj∆φ) hi(rj+2, θk,φNφ , tn) (2.97) hi(rj, θk,φNφ+1, tn) = −60r2j∆φ2 ∆z(25∆z + 12rj∆φ) fi(rj, θk, zNz , tn) + 60r2j∆φ 2 ∆z(25∆z + 12rj∆φ) fi(rj, θk, zNz−1, tn) + ￿ 15∆z 4(25∆z + 12rj∆φ) − 1 4 ￿ hi(rj, θk,φNφ−4, tn) − ￿ 80∆z 3(25∆z + 12rj∆φ) + 5 3 ￿ hi(rj, θk,φNφ−3, tn) + ￿ 90∆z 25∆z + 12rj∆φ − 5 ￿ hi(rj, θk,φNφ−2, tn) − ￿ 240∆z 25∆z + 12rj∆φ + 10 ￿ hi(rj, θk,φNφ−1, tn) + ￿ 2075∆z 12(25∆z + 12rj∆φ) − 65 12 ￿ hi(rj, θk,φNφ , tn) − 5rj∆φ 2∆z 2∆r(25∆z + 12rj∆φ) hi(rj−2, θk,φNφ , tn) + 20rj∆φ2∆z ∆r(25∆z + 12rj∆φ) hi(rj−1, θk,φNφ , tn) − 20rj∆φ 2∆z ∆r(25∆z + 12rj∆φ) hi(rj+1, θk,φNφ , tn) + 5rj∆φ2∆z 2∆r(25∆z + 12rj∆φ) hi(rj+2, θk,φNφ , tn) (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, differing by θ = π, at θ-index (k+ Nθ2 )mod Nθ. This allows demarkation: fi(r0, θk, zl) = fi(r1, θ(k+Nθ2 ) mod Nθ , zl), (2.105) fi(r−1, θk, zl) = fi(r2, θ(k+Nθ2 ) mod Nθ , zl) (2.106) 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θ2 ) mod Nθ ,φ1), (2.107) gi(rj, θk,φ−1) = gi(rj, θ(k+Nθ2 ) mod Nθ ,φ2), (2.108) gi(rj, θk,φ−2) = gi(rj, θ(k+Nθ2 ) mod Nθ ,φ3), (2.109) for each reaction species, i, within the bottom cap, and hi(rj, θk,φ0) = hi(rj, θ(k+Nθ2 ) mod Nθ ,φ1), (2.110) hi(rj, θk,φ−1) = hi(rj, θ(k+Nθ2 ) mod Nθ ,φ2), (2.111) hi(rj, θk,φ−2) = hi(rj, θ(k+Nθ2 ) mod Nθ ,φ3), (2.112) 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 differences for diffusion in rod-cell geometry contain coefficients that vary in orders of magnitude, based on spatial dimension and grid point location, arising from nonautonomous coefficients in rod cell diffusion equations, (2.3), (2.4), and (2.5). Large coefficient variation means that diffusion in rod-cell geometry evolves on different time scales. The explicit forward Euler method can solve diffu- sion equations in rod cell geometry, but only stably when using a very short time 28 step, as large coefficient variation results in instability under longer time steps. Im- plicit solution methods, such as the backward Euler method and Crank-Nicolson method, can numerically solve diffusion equations in rod cell geometry using longer time steps than the forward Euler method, but are computationally expensive. Nat- urally, an explicit numerical method that evolves only within pertinent time scales, on individual spatial dimensions, at individual grid points, would allow larger time- stepping than the forward Euler method, while maintaining more computational efficiency 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 dimen- sion finite difference coefficient, 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 differences in z evolve under the same variable time step as finite differences in φ for uniform time stepping throughout matching conditions. Within the rod-cell body: cr,i,j,k,l = Di ·max{1, r−1j }, (2.129) cθ,i,j,k,l = Di · r−2j , (2.130) cz,i,j,k,l = ￿ Di · r−2j for l = 1 or l = Nz; Di otherwise (2.131) Within both caps: cr,i,j,k,m = Di ·max{1, 2r−1j }, (2.132) cθ,i,j,k,m = Di · r−2j sin2 φm, (2.133) cφ,i,j,k,m = Di ·max ￿ cosφm r2j sinφm , r−2j ￿ . (2.134) The base 2 multi-time-step forward Euler method employs Nt different 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 differences, at each each grid point, in each spatial variable, evolve in only one time-step, chosen such that: χ∆τs(cψ,i,j,k,l) =  1 if 2Nt−s ≥ αi · cψ,i,j,k,l > 2Nt−s−1; 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 difference evolutions under longer or shorter time steps, chosen for each reaction species, i, to ensure stability and maximum computational efficiency. Time steps, in both caps, at each grid point, for each spatial variable, are chosen such that: χ∆τs(cψ,i,j,k,m) =  1 if 2Nt−s ≥ αi · cψ,i,j,k,m > 2Nt−s−1; 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 differences, 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)) = ∆τs ·Di ￿￿ 1− δl1 + δlNz 2 ￿ · ￿χ∆τs(Di ·max{1, r−1j }) ·∆2rfi(rj, θk, zl, tn) + χ∆τs(Di ·max{1, r−1j }) · 1 rj ∆1rfi(rj, θk, zl, tn) + χ∆τs(Di · r−2j ) · 1 r2j ∆2θfi(rj, θk, zl, tn) + (δl1 + δlNz) · χ∆τs(Di · r−2j ) ·∆2zfi(rj, θk, zl, tn) + (1− δl1 − δlNz) · χ∆τs(Di) ·∆2zfi(rj, θk, zl, tn) ￿￿ 30 +∆τs ·Di ￿￿ 1− δmNφ 2 ￿ · ￿χ∆τs(Di ·max{1, 2r−1j }) ·∆2rgi(rj, θk,φm, tn) + χ∆τs(Di ·max{1, 2r−1j }) · 2 rj ∆1rgi(rj, θk,φm, tn) + χ∆τs(Di · r−2j sin−2 φm) · 1 r2j sin 2 φm ∆2θgi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · cosφm r2j sinφm ∆1φgi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · 1 r2j ∆2φgi(rj, θk,φm, tn) ￿￿ +∆τs ·Di ￿￿ 1− δmNφ 2 ￿ · ￿χ∆τs(Di ·max{1, 2r−1j }) ·∆2rhi(rj, θk,φm, tn) + χ∆τs(Di ·max{1, 2r−1j }) · 2 rj ∆1rhi(rj, θk,φm, tn) + χ∆τs(Di · r−2j sin−2 φm) · 1 r2j sin 2 φm ∆2θhi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · cosφm r2j sinφm ∆1φhi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · 1 r2j ∆2φhi(rj, θk,φm, tn) ￿￿ (2.138) There are NrNθ(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,sfi,t = ∆sfi,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 ∆sfi,t is a column vector with Ngrid entries that correspond to the rod-cell geometry base 2 multi-time-step forward Euler finite difference 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 difference evolutions are calculated by initially incorporating only the finite difference evolutions that occur during the shortest time step, ∆t1, keeping constant all grid point values unaffected by finite difference 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 difference 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 difference evolutions occur. Thus, fi,t+∆t2 = ((I + Fi,1)(I + Fi,1) + Fi,2) fi,t. (2.141) Similarly, for all finite difference evolutions covering ∆t3: fi,t+∆t3 = ￿￿ (I + Fi,1) 2 + Fi,2 ￿ ￿ (I + Fi,1) 2 + Fi,2 ￿ + Fi,3 ￿ fi,t. (2.142) Continuing the iterative process, covering all finite difference evolutions under the base 2 multi-time-step forward Euler macro time step, ∆t: fi,t+∆t = ￿￿. . .￿￿(I + Fi,1)2 + Fi,2￿2 + Fi,3￿2 . . .￿2 + Fi,Nt−1 ￿2 + Fi,Nt  fi,t. (2.143) Ai, the complete base 2 multi-time-step forward Euler matrix, covering all finite difference evolutions over the macro time step ∆t, satisfies: fi,t+∆t = Aifi,t, (2.144) and is calculated: Ai = ￿￿. . .￿￿(I + Fi,1)2 + Fi,2￿2 + Fi,3￿2 . . .￿2 + Fi,Nt−1 ￿2 + Fi,Nt  , (2.145) where I is a Ngrid ×Ngrid sparse identity matrix. 2.4.4 Computational Efficiency of the Base 2 Multi-Time- Step Forward Euler Method Incorporating finite difference values only within pertinent time scales, while utiliz- ing 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 = ￿￿￿ . . . (((I + Fi,1) n + Fi,2) n + Fi,3) n . . . ￿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 differences within the next time interval, as op- posed 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 efficiency, 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, conver- gent 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−8s, 10−7s,. . . , 10−2s. Then, I calculated the l2 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µm2s−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−6s, 10−6s, 10−7s, and 10−3s. 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µm2s−1, Nglobal = 5, and αi = 1. Method time step BE CN FE MTS 10−8s 2.1 · 10−5 3.5 · 10−8 2.1 · 10−5 3.5 · 10−9 10−7s 1.9 · 10−3 3.4 · 10−5 2.2 · 10−3 3.4 · 10−7 10−6s .09 .02 .24 2.1 · 10−5 10−5s .34 - - 1.4 · 10−4 10−4s .68 - - 1.6 · 10−3 10−3s 1.6 - - .02 10−2s 5.0 - - .27 The base 2 multi-time-step forward Euler method maintains computational efficiency by evolving over a relatively long macro time step, while maintaing spar- sity. For example, when Di = 1µm2s−1, Nglobal = 5, and αi = 1, the base 2 multi-time-step forward Euler method’s maximum appropriate macro time step is 10−3s, while the forward Euler method’s maximum appropriate time step is 10−7s. Additionally, the complete base 2 multi-time-step forward Euler finite difference ma- trix, Ai, contains only about an order of magnitude more non-zero matrix elements than the analogous forward Euler finite difference 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 efficient, when cal- culating diffusion 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 difference 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 for- ward Euler and base 2 multi-time-step forward Euler finite difference 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. Calcu- lations are from rod-cells with Di = 1µm2s−1, under a grid with Nglobal = 5, αi = 1, and base 2 multi-time-step forward Euler macro time step ∆t = 10−3s. Rod-cell length 2µm 4µm 6µm 8µm 10µm Total Matrix Elements 4.4 · 106 16.8 · 106 37.2 · 106 65.6 · 106 102.0 · 106 FE NZME. 2.9 · 104 4.9 · 104 7.0 · 104 9.0 · 104 11.0 · 104 MTS NZME. 9.2 · 105 11.0 · 105 12.7 · 105 14.5 · 105 16.2 · 105 FE fraction NZME. 6.5 · 10−3 2.9 · 10−3 1.9 · 10−3 1.4 · 10−3 1.1 · 10−3 MTS fraction NZME. 21.0 · 10−2 6.5 · 10−2 3.4 · 10−2 2.2 · 10−2 1.6 · 10−2 The base 2 multi-time-step forward Euler method for solving diffusion equa- tions in rod-cell geometry is computationally more efficient than backward Euler or Crank-Nicolson implicit methods, as it implements sparse matrix multiplica- tion rather than more computationally expensive implicit matrix solving methods, while covering a larger time interval per solution calculation. For example, when Di = 1µm2s−1 and Nglobal = 5, the backward Euler and Crank-Nicolson methods’ maximum appropriate time step is 10−6s. Comparatively, under the same grid, the base 2 multi-time-step forward Euler method’s maximum appropriate macro time step is 10−3s, 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 different 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 diffusion equations by matrix multiplication is computationally expensive, but reduces solution com- putational 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µm2s−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µm2s−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 car- rying 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µm2s−1, Nglobal = 5, and αi = 1. Rod-cell length Method 2µm 4µm 6µm 8µm 10µm BE Implicit .6s .9s 1.3s 1.6s 2.1s BEM 9.9s 33.5s 139.9s 596.5s 864.7s CN Implicit .6s .9s 1.3s 1.6s 2.1s CNM 9.7s 35.0s 143.4s 614.7s 892.2s FE .6s .9s 1.3s 1.6s 2.1s MTS 5.7s 7.1s 8.2s 9.8s 11.1s Table 2.4: Computational time to calculate one numerical second using various methods. Definitions and parameter values are as in Table 2.3 . Rod-cell length Method 2µm 4µm 6µm 8µm 10µm BE Implicit 1.6 · 105s 2.5 · 105s 8.4 · 105s 1.6 · 106s 3.0 · 106s BEM 1.3 · 104s 5.0 · 104s 1.1 · 105s 3.3 · 105s 6.0 · 105s CN Implicit 1.6 · 105s 2.5 · 105s 9.1 · 105s 2.1 · 106s 2.9 · 106s CNM 1.3 · 104s 5.8 · 104s 1.1 · 105s 3.1 · 105s 5.9 · 105s FE 7.0 · 102s 1.2 · 103s 1.7 · 103s 3.1 · 103s 3.8 · 103s MTS 2.3s 3.3s 3.9s 4.4s 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 diffusion coefficients, 36 Di. Of the methods shown, the base 2 multi-time-step forward Euler method is clearly the most computationally efficient method for calculating diffusion in rod- cells 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 Diffusion in rod-cell geometry is a closed system. The total amount of each reaction species must be conserved during each time step of numerical diffusion. 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 vTfi,t. The concentration of reaction species i, at time t+∆t is Aifi,t, with total mass vTAifi,t. The base 2 multi-time-step forward Euler Method is conserving, for each reaction species i, if and only if: vTAifi,t = v Tfi,t, , (2.147) for any fi,t. vTAifi,t = v Tfi,t ￿ vTAifi,t − vTfi,t = 0 (2.148) ￿ vT(Aifi,t − fi,t) = 0 (2.149) ￿ vT(Ai − I)fi,t = 0 (2.150) ￿ fTi,t(Ai − I)Tv = 0 (2.151) 37 ￿ fTi,t(A T i − IT)v = 0 (2.152) ￿ fTi,t(A T i − I)v = 0. (2.153) If 1 is an eigenvalue of ATi with corresponding eigenvector, v, then (ATi − I)v = 0. (2.154) ￿ fTi,t(A T i − I)v = 0, (2.155) ￿ vTAifi,t = v Tfi,t for any fi,t. Thus, if ATi has eigenvalue, 1, and strictly positive or strictly negative (indi- cating a possible vector of concentrations) corresponding eigenvector, v, then the total amount of reaction species i is conserved during each time step of numerical diffusion, along v. The eigenvalue, 1, of ATi and corresponding strictly positive or strictly nega- tive eigenvalue, v, are verified using the MATLAB function “eigs”, which numeri- cally calculates eigenvalues of sparse matrices. If ATi has unique eigenvalue, 1, which is determined using the MATLAB func- tion “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￿Ngrid p=1 vp v, (2.156) 38 and satisfies: vTvolumeAifi,t = v T volumefi,t, (2.157) vTfi,t = total amount of reaction species, i, and (2.158) Ngrid￿ p=1 vvolume,p = rod-cell volume, (2.159) 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 diffusion 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 different grids, (Nglobal = {5, 6, 7, 8, 9, 10}). There are four grid points com- mon to all six grids: x1(r-index = Nr, θ-index = 1, z-index = 1), x2(r-index = Nr, θ-index = 1 + Nθ 2 , z-index = 1), x3(r-index = Nr, θ-index = 1, z-index = Nz), and x4(r-index = Nr, θ-index = 1+ Nθ 2 , z-index = Nz). A non-negative test function in cartesian coordinates, varying in all spatial directions, with differing 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, diffusion coefficient, Di = 1µm2s−1, and under the time step fine-tuning parameter, αi = 1. The ordered vector containing the numeric- volume of each grid point, under Nglobal grid points, is denoted as v Nglobal volume. For each Nglobal ∈ {5, 6, 7, 8, 9, 10}, vNglobalvolume · 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 differing total masses, when measuring spatial convergence. Numerical solution concentration differences, |f 10i,t (xj)− fni,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−4s. Numerical solutions differ at each grid point, x1, x2, x3 and x4, but numerical concentration differences are identical, up to numerical round- off, 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 differences for each Nglobal ∈ {5, 6, 7, 8, 9, 10}, at each numeric-volume, vNglobalvolume(xj), for grid points x1 and x4 and for grid points x2 and x3. Figure 2.1 also shows log plots of numerical concentration differences for each Nglobal ∈ {5, 6, 7, 8, 9}, at numeric-volume differ- ences, |v10volume(xj)− vnvolume(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 differences at volume differences shows the relationship: |f 10i,t (xj)− fni,t(xj)| = 17.9 · |v10volume(xj)− vnvolume(xj)|1.1, (2.161) at grid points x1 and x4, and |f 10i,t (xj)− fni,t(xj)| = 7.7 · |v10volume(xj)− vnvolume(xj)|0.9, (2.162) at grid points x2 and x3, n ∈ {5, 6, 7, 8, 9}. Therefore, concentration differences converge locally in numeric-volume differences at approximately order 1.0. Numeric- volume differences are approximately the difference of cubed grid spacing, hence concentration differences converge locally in grid spacing at approximately order 3, above the minimum spatial convergence order of 2, as predicted by finite difference formula derivations. 40 2 4 6 8 10 12 x 10−4 0 0.002 0.004 0.006 0.008 0.01 v volume n (xj) (µm 3) |f i,t10 (x j )−f i,tn (x j )| ( µm − 3 ) (a) −4.2 −4 −3.8 −3.6 −3.4 −3.2 −3 −3.2 −3 −2.8 −2.6 −2.4 −2.2 −2 log10|vvolume10 (xj)−vvolumen (xj)| lo g 1 0|f i ,t10 (x j )−f i,tn (x j )| (b) 2 4 6 8 10 12 x 10−4 0 0.002 0.004 0.006 0.008 0.01 0.012 v volume n (xj) (µm 3) |f i,t10 (x j )−f i,tn (x j )| ( µm − 3 ) (c) −4.2 −4 −3.8 −3.6 −3.4 −3.2 −3 −3 −2.8 −2.6 −2.4 −2.2 −2 log10|vvolume10 (xj)−vvolumen (xj)| lo g 1 0|f i ,t10 (x j )−f i,tn (x j )| (d) Figure 2.1: Spatial convergence of the base 2 multi-time-step forward Euler method in rod-cell geometry. Numerical concentration differences, 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 differences, for each Nglobal ∈ {5, 6, 7, 8, 9}, at logs of grid point volume differences, 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 For- ward Euler Method Numerical solutions to the system of diffusion 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 time step, ∆t, and time step fine-tuning parameter, αi = 1, is denoted as f∆ti,t . To test for temporal convergence, numerical solutions, f∆ti,t , were found over the time interval, 10−3s, using macro time steps, ∆tn = 2−n ·10−3s, 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 diffusion coefficient, Di = 1µm2s−1. Differences in numerical solutions from the numerical solution found under macro time step, ∆t5, ||f∆t5i,t − f∆tni,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 differences in numerical solutions from the numerical solution found under macro time step, ∆t5, at the log of differences in macro time step from the macro time step, ∆t5, shown in Figure 2.2, shows the relationship: ||f∆t5i,t − f∆tni,t || = 2.5 · 103 · |∆t5 −∆tn|1.0, (2.163) for all n ∈ {0, 1, 2, 3, 4}. Hence, numerical solutions converge temporally at approx- imately O(∆t), as expected for a forward Euler method. 42 −15 −14 −13 −12 −11 −10 −9 0 0.5 1 1.5 2 2.5 log2(∆ t=2 −n ⋅10−3) ||f i ,t∆ t= 2− 5 ⋅ 10 − 3 − f i, t∆ t= 2− n ⋅ 10 − 3 || (a) −15 −14 −13 −12 −11 −10 −5 −4 −3 −2 −1 0 1 2 log2|∆ t=2−5⋅10−3−∆ t=2−5⋅10−n| lo g 2 ||f i ,t∆ t= 2− 5 ⋅ 10 − 3 − f i, t∆ t= 2− n ⋅ 10 − 3 || (b) Figure 2.2: Temporal convergence of the base 2 multi-time-step forward Euler method. Differences in numerical solutions found under all macro time steps, ∆tn = 2−n · 10−3s, 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 differences 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 differences 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 diffusion equation in three dimensions, on an infinite domain, starting with δ function at (x = x0, y = y0, z = z0, t = 0) is: u(x, y, z, t) = 1 (4πDt)3/2 exp ￿ −(x− x0) 2 + (y − y0)2 + (z − z0)2 4Dt ￿ . (2.164) Numerical convergence studies demonstrate both spatial, in Section 2.4.6, and tem- poral, in Section 2.4.7, convergence of numerical solutions to the system of diffusion equations, Equation 2.2, in rod-cell geometry, using the base 2 multi-time-step for- ward Euler method. Beyond spatial and temporal convergence, numerical diffusion 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 sufficiently far from the rod-cell membrane, with sufficiently small variance, such that numerical solutions at rod-cell boundary grid points have sufficiently small values, for some sufficiently small time interval. In this regime, numerical solutions on the bounded rod-cell domain should closely approximate the infinite-domain so- lution. To test convergence of numerical solutions to u(x, y, z, t), a numerical so- lution was found using the sufficiently small macro time step, ∆t = 10−4s, over the sufficiently small time interval, [0, 5 · 10−3s], starting from the initial condition u(x, y, z, t0 = 4 · 10−3s), centered at the interface of the rod-cell bottom cap and body, sufficiently 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 diffusion coefficient, D = Di = 1µm2s−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−3s. Relative differences 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 for- ward Euler diffusion 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 sufficiently consistent conditions. Addition- ally, 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). 0 1 2 3 4 5 x 10−3 4.5 5 5.5 6 6.5 7 7.5 8 8.5 x 10−3 n⋅∆ t (s) σ (t) (a) 0 1 2 3 4 5 x 10−3 0 1 2 3 4 5 6 x 10 −3 n⋅∆ t (s) ||u t 0 + n ⋅ ∆ t− f i, n⋅ ∆ t||⋅| |u t 0+ n ⋅ ∆ t||− 1 (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 time- dependent 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−4s. 2.4.9 Diffusion 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 diffusion require numerical diffusion solutions that are restricted to the rod-cell boundary. Within the restriction of the system of diffusion equations, Eq. 2.2, to the rod cell boundary, there is no diffusion in the r direction, and diffusion in other spatial variables only occurs at the rod cell boundary, where r = R. Surface diffusion is described by: ∂f̄i(R, θ, z, t) ∂t = Di ￿ 1 R2 ∂2f̄i(R, θ, z, t) ∂θ2 + ∂2f̄i(R, θ, z, t) ∂z2 ￿ , (2.167) on the rod-cell body boundary, ∂ḡi(R, θ,φ, t) ∂t =Di ￿ 1 R2 sin2 φ ∂2ḡi(R, θ,φ, t) ∂θ2 + cosφ R2 sinφ ∂ḡi(R, θ,φ, t) ∂φ + 1 R2 ∂2ḡi(R, θ,φ, t) ∂φ2 ￿ , (2.168) 45 on the bottom cap boundary, and ∂h̄i(R, θ,φ, t) ∂t =Di ￿ 1 R2 sin2 φ ∂2h̄i(R, θ,φ, t) ∂θ2 + cosφ R2 sinφ ∂h̄i(R, θ,φ, t) ∂φ + 1 R2 ∂2h̄i(R, θ,φ, t) ∂φ2 ￿ , (2.169) on the top cap boundary, for each membrane-bound reaction species i. Boundary conditions and cap-body interface matching conditions for diffusion on the rod-cell boundary are a restriction, to r = R, of boundary conditions for diffusion in rod-cell geometry. Numerical diffusion solutions on the rod-cell boundary are found using the rod-cell grid, as defined in Section 2.4, pertinent rod-cell finite difference 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-time- step forward Euler method, as described in Section 2.4.3. Thus, base 2 multi-time- step forward Euler finite differences for diffusion on the rod-cell boundary:￿ f̄i(rj, θk, zl, tn +∆τs) + ḡi(rj, θk,φm, tn +∆τs) + h̄i(rj, θk,φm, tn +∆τs) ￿ − ￿f̄i(rj, θk, zl, tn) + ḡi(rj, θk,φm, tn) + h̄i(rj, θk,φm, tn)￿ = ∆τs ·Di ￿￿ 1− δl1 + δlNz 2 ￿ · δjNr · ￿ χ∆τs(Di · r−2j ) · 1 r2j ∆2θf̄i(rj, θk, zl, tn) + (δl1 + δlNz) · χ∆τs(Di · r−2j ) ·∆2zf̄i(rj, θk, zl, tn) + (1− δl1 − δlNz) · χ∆τs(Di) ·∆2zf̄i(rj, θk, zl, tn) ￿￿ +∆τs ·Di ￿￿ 1− δmNφ 2 ￿ · δjNr · ￿ χ∆τs(Di · r−2j sin−2 φm) · 1 r2j sin 2 φm ∆2θḡi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · cosφm r2j sinφm ∆1φḡi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · 1 r2j ∆2φḡi(rj, θk,φm, tn) ￿￿ 46 +∆τs ·Di ￿￿ 1− δmNφ 2 ￿ · δjNr · ￿ χ∆τs(Di · r−2j sin−2 φm) · 1 r2j sin 2 φm ∆2θh̄i(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · cosφm r2j sinφm ∆1φh̄i(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · 1 r2j ∆2φh̄i(rj, θk,φm, tn) ￿￿ (2.170) F̄i,s is the Ngrid ×Ngrid sparse matrix, such that: F̄i,sf̄i,t = ∆sf̄i,t, (2.171) where f̄i,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 ∆sf̄i,t is a column vector with Ngrid entries that correspond to the rod-cell boundary base 2 multi-time-step forward Euler finite difference value:￿ f̄i(rj, θk, zl, tn +∆τs) + ḡi(rj, θk,φm, tn +∆τs) + h̄i(rj, θk,φm, tn +∆τs) ￿ − ￿f̄i(rj, θk, zl, tn) + ḡi(rj, θk,φm, tn) + h̄i(rj, θk,φm, tn)￿ , for the time step, ∆τs, of membrane-bound reaction species, i, at time, t, at each grid point, ordered from 1 to Ngrid. Āi is the complete rod-cell boundary base 2 multi-time-step forward Euler matrix, covering all finite difference evolutions over the macro time step, ∆t. Āi satisfies: f̄i,t+∆t = Āif̄i,t, (2.172) and is calculated: Āi = ￿. . .￿￿￿I + F̄i,1￿2 + F̄i,2￿2 + F̄i,3￿2 . . .￿2 + F̄i,Nt−1 2 + F̄i,Nt  , (2.173) where I is a Ngrid ×Ngrid sparse identity matrix. If Ăi = Āi, excluding all dimensions that do not correspond to grid points on the rod-cell boundary, then as shown for rod-cell diffusion in Section 2.4.5, conserva- tion during each time step of numerical diffusion, on the boundary, is determined by eigenvalue, 1, with corresponding strictly positive or strictly negative eigenvector, 47 v̆, of ĂTi . If 1 is a unique eigenvalue of Ă T i , then an ordered vector, of which each 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￿Nθ(2Nφ+Nz−2) p=1 vp v̆, (2.174) where v̆p is the pth element of v̆. The MATLAB function “eigs” was used to verify unique eigenvalue, 1, and corre- sponding strictly positive or strictly negative eigenvector, v̆, of ĂTi . Numerical solutions to the system of diffusion equations, Equation 2.2, in rod- cell 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. Diffusion on the rod-cell boundary is a restriction of rod-cell diffusion, lacking only diffusion in the r direc- tion, 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 diffusion equations, Equation 2.1, are numerically solved at each grid point in rod-cell geometry us- ing 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 diffusion solutions to form a complete solution to the reaction-diffusion system, Equation 2.1. This step splitting procedure is reasonable if reactions are relatively slow compared to the diffusive time scale. Thus, macro time steps are chosen to be sufficiently short, such that concentration changes from reactions are not significantly greater than concentration changes from diffusion, over the macro time step. 48 Chapter 3 Numerical Solutions to Reaction-Diffusion 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 diffusion equations in growing and dividing rod-cell geometries. 3.1 Previously Published Results and Numerical Methods for Solving Reaction-Diffusion Sys- tems 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-diffusion sys- tem 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 di- mensional growing domain. Kulesa et al. demonstrated the spatial patterning of alligator Teeth Primordia, through a one dimensional reaction-diffusion 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 split- ting in a growing one dimensional domain (Meinhardt and de Boer, 2001). Nu- merically, 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 possi- ble mechanism for Min protein distribution in daughter cells (Ventura and Sourjik, 2011). Tostevin and Howard’s stochastic, one dimensional, Min system model ac- complished dynamic septation by sequentially reducing the probability of diffusion across the domain mid point (Tostevin and Howard, 2006). Sengupta and Ruten- berg deterministically modeled dynamic partitioning between daughter cells, using the Huang et al model, in a three dimensional cylindrical representation of Escheri- cia Coli, through sequential closure of diffusion 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 numeri- cal method for rod-cell division, similar to the Sengupta and Rutenberg method, though permitting diffusion of membrane-bound reaction species, with variable sep- tation locations. 3.2 Numerically Solving Reaction-Diffusion Equa- tions 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-diffusion 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 = 2LNNglobal + 1, is set, and will be used for all rod-cell body lengths during growth. Next, reaction-diffusion 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: ∆z0 = L0 Nz − 1 , (3.1) ∆z1 = L1 Nz − 1 = L0 +∆L1 Nz − 1 . (3.2) The value of each reaction species, i, is unchanged at every grid point, but the rod- cell has grown in length. Reaction-diffusion equations are then numerically solved, using the base 2 multi-time-step forward Euler and RK4 methods, under the rod- cell body length L1. The process is repeated until the final rod-cell body length, LN , has been attained. Then, reaction-diffusion 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 − L0 = N￿ p=1 ∆Lp (3.3) 3.2.3 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 efficient 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, re- sulting in smoother growth. Each growth method spatially distorts solutions. Growth by adding cylin- drical 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 relatively small, of O( ∆LpNz−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 efficient as growth by adding cylin- drical 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 cylindri- cal 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-diffusion 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 diffusion 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￿ p=1 ∆Lqp. (3.4) 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 point within the rod-cell body, is of O( ∆L q p 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 dupli- cated 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 diffusion barrier for bulk reaction species and as a domain extension for membrane 53 bound reaction species. Initially, there is no invagination. Then, at each invagi- nation step, q ∈ {1, 2, . . . , Nr}, the boundary invaginates to R − (2(q−1)+1)∆r2 , the 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 be- tween 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 Diffusion Equations in Dividing Rod-Cells Solutions to the diffusion parts of reaction-diffusion equations, in dividing rod-cells, are found using the base 2 multi-time-step forward Euler method, described in Chap- ter 2. All spatial finite difference formulas, from Section 2.4.1, and all boundary conditions, from Section 2.4.2, for solving diffusion equations in rod-cell geometry are applicable for solving diffusion equations in dividing rod-cells, but extra condi- tions 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 difference 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 diffusion is open in z across zd and zd+1, while maintaining the proper magnitude of diffusivity, the total diffusion in r and θ, at zd and zd+1, is half the diffusion in r and θ at zd, and half the diffusion in r and θ at zd+1. At each invagination step, q, the boundary has invaginated in the r direc- tion past grid points of r-index Nr − (q − 1). Thus, at grid points of r-index ≥Nr − (q− 1), diffusion is closed in the positive z direction at zd, and closed in the negative z-direction at zd+1. 54 Second order finite differences in z, at zd and zd+1 are: ∂2f(rj, θk, zd, tn) ∂z2 ≈f(rj, θk, zd−1, t) ∆z2 − 2f(rj, θk, zd, t) ∆z2 + f(rj, θk, zd+2, t) ∆z2 , (3.5) ∂2f(rj, θk, zd+1, tn) ∂z2 ≈f(rj, θk, zd−1, t) ∆z2 − 2f(rj, θk, zd+1, t) ∆z2 + f(rj, θk, zd+2, t) ∆z2 . (3.6) 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 diffusion equations in dividing rod- cells, in conjunction with the base 2 multi-time-step forward Euler method for solv- ing diffusion equations in rod-cell geometry, from Section 2.4.3, result in base 2 multi-time-step forward Euler finite differences for solving diffusion 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)) = ∆τs ·Di ￿￿ 1− δl1 − δlNz 2 − δld − δl(d+1) ￿ ·￿ χ∆τs(Di ·max{1, r−1j }) ·∆2rfi(rj, θk, zl, tn) + χ∆τs(Di ·max{1, r−1j }) · 1 rj ∆1rfi(rj, θk, zl, tn) + χ∆τs(Di · r−2j ) · 1 r2j ∆2θfi(rj, θk, zl, tn) + (δl1 + δlNz) · χ∆τs(Di · r−2j ) ·∆2zfi(rj, θk, zl, tn) + (1− δl1 − δlNz) · χ∆τs(Di) ·∆2zfi(rj, θk, zl, tn) ￿ 55 + δld 2 · (1 +H[j − (Nr − (q − 1))]) · ￿ χ∆τs(Di ·max{1, r−1j }) ·∆2rfi(rj, θk, zl, tn) + χ∆τs(Di ·max{1, r−1j }) · 1 rj ∆1rfi(rj, θk, zl, tn) + χ∆τs(Di · r−2j ) · 1 r2j ∆2θfi(rj, θk, zl, tn) ￿ + δld 2 · (1−H[j − (Nr − (q − 1))]) · ￿ χ∆τs(Di ·max{1, r−1j }) ·∆2rfi(rj, θk, zl+1, tn) + χ∆τs(Di ·max{1, r−1j }) · 1 rj ∆1rfi(rj, θk, zl+1, tn) + χ∆τs(Di · r−2j ) · 1 r2j ∆2θfi(rj, θk, zl+1, tn) ￿ + δl(d+1) 2 · (1 +H[j − (Nr − (q − 1))]) · ￿ χ∆τs(Di ·max{1, r−1j }) ·∆2rfi(rj, θk, zl, tn) + χ∆τs(Di ·max{1, r−1j }) · 1 rj ∆1rfi(rj, θk, zl, tn) + χ∆τs(Di · r−2j ) · 1 r2j ∆2θfi(rj, θk, zl, tn) ￿ + δl(d+1) 2 · (1−H[j − (Nr − (q − 1))]) · ￿ χ∆τs(Di ·max{1, r−1j }) ·∆2rfi(rj, θk, zl−1, tn) + χ∆τs(Di ·max{1, r−1j }) · 1 rj ∆1rfi(rj, θk, zl−1, tn) + χ∆τs(Di · r−2j ) · 1 r2j ∆2θfi(rj, θk, zl−1, tn) ￿ + δld · χ∆τs(Di) · ￿ fi(rj, θk, zl−1, tn) ∆z2 − 2fi(rj, θk, zl, tn) ∆z2 + fi(rj, θk, zl+2, tn) ∆z2 ￿ + δl(d+1) · χ∆τs(Di) · ￿ fi(rj, θk, zl−2, tn) ∆z2 − 2fi(rj, θk, zl, tn) ∆z2 + fi(rj, θk, zl+1, tn) ∆z2 ￿￿ 56 +∆τs ·Di ￿￿ 1− δmNφ 2 ￿ · ￿χ∆τs(Di ·max{1, 2r−1j }) ·∆2rgi(rj, θk,φm, tn) + χ∆τs(Di ·max{1, 2r−1j }) · 2 rj ∆1rgi(rj, θk,φm, tn) + χ∆τs(Di · r−2j sin−2 φm) · 1 r2j sin 2 φm ∆2θgi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · cosφm r2j sinφm ∆1φgi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · 1 r2j ∆2φgi(rj, θk,φm, tn) ￿￿ +∆τs ·Di ￿￿ 1− δmNφ 2 ￿ · ￿χ∆τs(Di ·max{1, 2r−1j }) ·∆2rhi(rj, θk,φm, tn) + χ∆τs(Di ·max{1, 2r−1j }) · 2 rj ∆1rhi(rj, θk,φm, tn) + χ∆τs(Di · r−2j sin−2 φm) · 1 r2j sin 2 φm ∆2θhi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · cosφm r2j sinφm ∆1φhi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · 1 r2j ∆2φhi(rj, θk,φm, tn) ￿￿ (3.9) where H is the discrete Heaviside function. Using dividing rod-cell base 2 multi-time-step forward Euler finite differences, Ai, the complete base 2 multi-time-step forward Euler matrix, covering all finite difference evolutions, over the macro time step, ∆t, is calculated as in Section 2.4.3. Conservation during each time step of numerical diffusion, in the dividing rod-cell, is determined using the method of Section 2.4.5. Complete solutions to reaction- diffusion 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 Diffusion Equations on Dividing Rod-Cell Boundaries The dividing rod-cell boundary is essentially the same as the the rod-cell bound- ary, with added invagination in the r direction at zd and zd+1. Solving diffusion equations on the dividing rod-cell boundary follows directly from solving diffusion equations on the rod-cell boundary, as in Section 2.4.9, with added conditions at invaginated boundary points. Membrane bound proteins can diffuse around the leading edge of the invagi- nated 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 r − (q−2)rN r − (q−1)rN r − (q−1)rN r − (q−2) zd zd+1 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 im- plementation of symmetric five point finite difference formulas in r, shown in Section 58 2.4.1. To allow diffusion around the leading edge of the invaginated boundary, con- ditions: f̄i(rNr−q, θk, zd) = f̄i(rNr−(q−1), θk, zd+1), (3.10) f̄i(rNr−q−1, θk, zd) = f̄i(rNr−(q−2), θk, zd+1), (3.11) f̄i(rNr−q, θk, zd+1) = f̄i(rNr−(q−1), θk, zd), (3.12) f̄i(rNr−q−1, θk, zd+1) = f̄i(rNr−(q−2), θk, zd), (3.13) are set in finite difference formulas in r, for q ∈ {1, 2, . . . , Nr−1}. When q = Nr, dif- fusion is closed across the invaginated boundary, hence finite difference conditions, in r, follow those for bulk rod-cell diffusion, 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 direc- tion across zd+1, hence second order finite difference 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 for- ward Euler finite differences for solving diffusion equations on the rod-cell boundary, shown in 2.4.9, result in base 2 multi-time-step forward Euler finite differences for solving diffusion equations on the dividing rod cell boundary:￿ f̄i(rj, θk, zl, tn +∆τs) + ḡi(rj, θk,φm, tn +∆τs) + h̄i(rj, θk,φm, tn +∆τs) ￿ − ￿f̄i(rj, θk, zl, tn) + ḡi(rj, θk,φm, tn) + h̄i(rj, θk,φm, tn)￿ = ∆τs ·Di ￿￿ 1− δl1 + δlNz 2 − δld − δl(d+1) ￿ · δjNr · ￿ χ∆τs(Di · r−2j ) · 1 r2j ∆2θf̄i(rj, θk, zl, tn) + (δl1 + δlNz) · χ∆τs(Di · r−2j ) ·∆2zf̄i(rj, θk, zl, tn) + (1− δl1 − δlNz) · χ∆τs(Di) ·∆2zf̄i(rj, θk, zl, tn) ￿ 59 + ￿ δld + δl(d+1) ￿ ·H[j − (Nr − (q − 1))]·￿ χ∆τs(Di ·max{1, r−1j }) ·∆2r f̄i(rj, θk, zl, tn) + χ∆τs(Di ·max{1, r−1j }) · 1 rj ∆1r f̄i(rj, θk, zl, tn) ￿ (3.14) + δjNr · δld · χ∆τs(Di) · ￿ f̄i(rj, θk, zl−1, tn) ∆z2 − 2f̄i(rj, θk, zl, tn) ∆z2 + f̄i(rj, θk, zl+2, tn) ∆z2 ￿ + δjNr · δl(d+1) · χ∆τs(Di) · ￿ f̄i(rj, θk, zl−2, tn) ∆z2 − 2f̄i(rj, θk, zl, tn) ∆z2 + f̄i(rj, θk, zl+1, tn) ∆z2 ￿￿ +∆τs ·Di ￿￿ 1− δmNφ 2 ￿ · δjNr · ￿ χ∆τs(Di · r−2j sin−2 φm) · 1 r2j sin 2 φm ∆2θḡi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · cosφm r2j sinφm ∆1φḡi(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · 1 r2j ∆2φḡi(rj, θk,φm, tn) ￿￿ +∆τs ·Di ￿￿ 1− δmNφ 2 ￿ · δjNr · ￿ χ∆τs(Di · r−2j sin−2 φm) · 1 r2j sin 2 φm ∆2θh̄i(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · cosφm r2j sinφm ∆1φh̄i(rj, θk,φm, tn) + χ∆τs(Di ·max{ cosφm r2j sinφm , r−2j }) · 1 r2j ∆2φh̄i(rj, θk,φm, tn) ￿￿ (3.15) Using dividing rod-cell boundary base 2 multi-time-step forward Euler finite differences, Āi, the complete base 2 multi-time-step forward Euler matrix, covering all finite difference evolutions, over the macro time step, ∆t, is calculated as in Section 2.4.9. Conservation during each time step of numerical diffusion, on the dividing rod-cell boundary, is determined using the method of Section 2.4.9. Com- plete solutions to reaction-diffusion 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-Diffusion Solutions in Rod-Cell Geometry To demonstrate the methods, developed in Chapters 2 and 3, for numerically solving reaction-diffusion 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 con- tinuously grows, are shown. Any reaction-diffusion model may be solved separately under static cell size, growth, or division to develop a more in depth understand- ing 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 cy- tosolic MinD:ATP complexes, ρD:ATP , the concentration of cytosolic MinD:ATP complexes, ρE, the concentration of cytosolic MinE, ρd, the concentration of mem- brane 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 and bind an ATP, forming MinD:ATP complexes, at the rate, σADP→ATPD . 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, mem- brane 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 diffusion on the boundary. Including diffusion on the boundary, the Huang et. al model becomes: ∂ρD:ADP ∂t =DD∇2ρD:ADP − σADP→ATPD ρD:ADP + δ(r −R)σdeρde, (4.1) ∂ρD:ATP ∂t =DD∇2ρD:ATP + σADP→ATPD ρD:ADP − δ(r −R)[σD + σdD(ρd + ρde)]ρD:ATP , (4.2) ∂ρE ∂t =DE∇2ρE + δ(r −R)σdeρde − δ(r −R)σEρdρE, (4.3) ∂ρd ∂t =Dd∇2 |R ρd − σEρdρE(R) + [σD + σdD(ρd + ρde)]ρD:ATP (R), (4.4) ∂ρde ∂t =Dde∇2 |R ρde − σdeρde + σEρdρE(R). (4.5) Parameters within the Huang et al. model have values: DD = DE =2.5µm 2s−1, (4.6) Dd = Dde =0µm 2s−1, (4.7) σADP→ATPD =1s −1, (4.8) σD =0.025µms −1, (4.9) σdD =0.0015µm 3s−1, (4.10) σde =0.7s −1, (4.11) σE =0.093µm 3s−1. (4.12) 63 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 uni- form 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 nor- mal 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 pa- rameter, 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 correc- tions were required, analogous to constant production of ρD:ADP . When negative MinD corrections were required, ρD:ADP and ρD:ATP 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. One- dimensional mean solution concentrations are shown for each reaction species: (a), ρD:ADP , (b), ρD:ATP , (c), ρD:ATP , 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 behav- iors, 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 diffusion coefficients are less than 2.5µm2s−1, when starting from the same initial condition (results not shown), indicating that larger cytosolic protein diffusivity may connect reaction re- gions sufficiently 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 concen- tration 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), (Cytryn- baum 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 solu- tions, (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 differing daughter cell behavior, a possible cause of the Fischer-Friedrich findings. No oscillatory behavior was observed when incorporat- ing diffusion on the rod-cell boundary, when boundary diffusion coefficients were chosen an order of magnitude smaller than cytosolic diffusion coefficients (results not shown). Membrane bound proteins likely diffuse in vivo. A future investigation of how and to what extent diffusion 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, rep- resenting 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, correc- tions 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:ATP , (c), ρD:ATP , 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 periodici- 68 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 efficient base 2 multi-time- step forward Euler method for solving reaction-diffusion 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, diffusivity, 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 effect of cell growth on bifurcation processes; and how divided rod- cell sub-domain behaviors differ, 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 mathemat- ically interesting Min system models, provide insight into new models, and produce biologically pertinent predictions. 69 Bibliography E. J. Crampin, E. A. Gaffney, and P. K. Maini. Reaction and diffusion on grow- ing 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é, and K. Karsten. Intra- and 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-diffusion 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. Mur- ray. 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 difference discretization for Poisson equation on a disk. Numer. Meth. Partial Diff. 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 mem- brane 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ón, 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 phos- pholipid vesicles regulated by ATP and MinE. PNAS, 99(10):6761–6766, may 2002. 72