Describing the movement of molecules in reduced-dimension models
1 Describing the movement of molecules in reduced-dimension models.
Natasha S. Savage. University of Liverpool, Liverpool, UK.
Email: [email protected]
Abstract
When addressing spatial biological questions using mathematical models, symmetries within biological systems are often exploited to simplify the problem by reducing its physical dimension. In reduced-dimension models molecular movement is restricted to the reduced dimension, changing the nature of molecular movement, leading to quantitatively and even qualitatively inaccurate results. Within this manuscript we discuss the condition under which restricted molecular movement in reduced-dimension models accurately approximate molecular movement in the full system. The condition depends on symmetries in geometry and pattern. For those systems which do not meet the condition, we present a general method for approximating unrestricted molecular movement in reduced-dimension models. We go on to use the ‘general method’ to derive a novel, mathematical robust, finite difference method for solving the 2D diffusion equation within a 1D reduced-dimension model.
Introduction
Simple models of complex systems are an invaluable tool for gaining conceptual insight into biological mechanism (1, 2). In spatial models a powerful simplification technique often used is to exploit symmetries within a biological system’s geometry and patterning to reduce the physical dimension of the problem. For example consider a single cell and the formation of a concentrated patch of membrane associated proteins, a polarity patch, this system has radial symmetry and so mechanisms controlling patch formation could be explored by considering a one dimensional (1D) slice through the centre of the patch, rather than considering the entire two dimensional (2D) membrane (3-8), Fig. 1A. When addressing questions which include cytoplasmic gradients, for example, a system’s dimension is often reduced from three to 2D (1, 9-13), or even 1D (1, 2, 14-16). An analogy for reduced-dimension models is the focal plane in microscopy: Analysis is performed on data acquired from a representative slice through the system, then results are inferred onto the entire cell or tissue. All spatial models contain a description of molecular movement. Molecular movement in reduced-dimension models is restricted to the focal plane, compare Fig. 1B-C. Ideally, molecular movement in reduced-dimension models would approximate molecules moving through the ‘focal plane’, Fig. 1B, we will show that this is only the case for a minority of reduced-dimension models. Without loss of generality, within this manuscript, we will consider the full system to be a 2D surface, a membrane, and the reduced-dimension model to be a 1D ring (Fig.1A). Molecular movement will be diffusive. Within this setting we will: 1) Derive the limited conditions under which the 1D diffusion equation accurately describes diffusive movement through the focal plane. 2) Show that, in general, the 1D diffusion equation produces quantitatively and even qualitatively inaccurate results when used in a reduced-dimension model. 3) Exploit the very symmetry used for dimension reduction to present a method for the approximation of molecular flux through the focal plane. We will then use the method to solve the 2D diffusion equation within the 1D reduced-dimension model, providing a 2 considerable extension to the set of accurate 1D reduced-dimension models. The scenarios to be used for illustrative purposes are; diffusion from a concentrated patch to homogeneity, estimating diffusion coefficients using FRAP analysis, and reaction-diffusion dynamics.
Results
To understand the inbuilt assumption of reduced-dimension models cover the 2D membrane with a regular mesh, Fig. 1D. The mesh-points along the 𝑥 -axis are labelled 𝑖 = 1,2, … , 𝑁 and are distance ∆𝑥 𝜇𝑚 apart, the mesh-points along the 𝑦 -axis are labelled 𝑗 = 1,2, … , 𝑀 and are distance ∆𝑦 𝜇𝑚 apart. Let 𝑢 𝑖,𝑗𝜏 represent the concentration of molecule 𝑢 on the membrane at point (𝑖, 𝑗) and time 𝜏 . The change in concentration of 𝑢 𝑖,𝑗𝜏 over time, as a result of diffusive movement, is described by the solution to the 2D diffusion equation. The explicit 2D finite difference scheme (2D-FDM) used to solve the 2D diffusion equation on a regular mesh, with diffusion coefficient 𝐷 𝜇𝑚 𝑠 −1 , is (17), 𝑢 𝑖,𝑗𝜏+1 = 𝑢 𝑖,𝑗𝜏 + ∆𝑡∆𝑥 𝐷(𝑢 𝑖−1,𝑗𝜏 − 2𝑢 𝑖,𝑗𝜏 + 𝑢 𝑖+1,𝑗𝜏 ) + ∆𝑡∆𝑦 𝐷(𝑢 𝑖,𝑗−1𝜏 − 2𝑢 𝑖,𝑗𝜏 + 𝑢 𝑖,𝑗+1𝜏 ) Note that to find 𝑢 𝑖,𝑗𝜏+1 the 2D-FDM calculates molecular movement along the 𝑥 -axis and 𝑦 -axis separately. To build a 1D reduced-dimension model describing the diffusive movement of 𝑢 on a 2D surface, the focal plane is set to run through the axis of symmetry of 𝑢 and the membrane, Fig. 1A. Molecule movement through the focal plane is then approximated by the solution to the 1D diffusion equation. Recall the 2D mesh, assume that the focal plane is set along the 𝑥 -axis at row 𝑗 = 𝐽 , Fig. 1D. The explicit 1D finite difference scheme (1D-FDM) used to solve the 1D diffusion equation is (17), 𝑢 𝑖,𝐽𝜏+1 = 𝑢 𝑖,𝐽𝜏 + ∆𝑡∆𝑥 𝐷(𝑢 𝑖−1,𝐽𝜏 − 2𝑢 𝑖,𝐽𝜏 + 𝑢 𝑖+1,𝐽𝜏 ) As the 1D reduced-dimension model only contains information about concentrations on row 𝐽 the 1D-FDM contains no terms for calculating the movement of 𝑢 along the 𝑦 -axis (compare the 2D and 1D-FDMs). Thus, molecular movement in the 1D reduced-dimension model is modelled as though it is restricted to the focal plane Fig. 1C. An inbuilt assumption of the 1D reduced-dimension model is that (𝑢 𝑖,𝐽−1𝜏 − 2𝑢 𝑖,𝐽𝜏 + 𝑢 𝑖,𝐽+1𝜏 ) = 0 , the number of molecules leaving the focal plane, −2𝑢 𝑖,𝐽𝜏 , is equal to the number of molecules entering it, 𝑢 𝑖,𝐽−1𝜏 + 𝑢 𝑖,𝐽+1𝜏 , at all points in space, 𝑖 , for all time, 𝜏 : There is zero flux through the focal plane. The zero flux assumption is only valid for a small subset of reduced-dimension models, those for which a mesh can be drawn such that the zero flux assumption holds. For example, consider the formation of a polarity patch on a spherical cell, one could construct a spherical mesh with a pole located at the centre of the patch, Fig. 1A. Because of the placing of the mesh concentrations 𝑢 𝑖,𝐽−1𝜏 , 𝑢 𝑖,𝐽𝜏 , and 𝑢 𝑖,𝐽+1𝜏 are equal and the zero flux assumption holds, Fig. 1E. Thus, molecular movement in a 1D reduced-dimension model, with focal plane at 𝑗 = 𝐽 , would be accurately described by the 1D diffusion equation. For the majority of cases the zero flux assumption does not hold. For example, the formation of a polarity patch along the body of an elongated cell, Fig. 1A. To solve the diffusion equation here a square mesh is constructed on the cell surface. As a polarity patch is radial and the mesh is square a focal plane cannot be found such that a 1D reduced-dimension model would obey the zero flux assumption, Fig. 1D: The standard 1D diffusion equation gives an inaccurate approximation of 2D molecular movement through the focal plane. 3 We now present a methodology that can be used to retain the significant reduction in model complexity afforded by dimension reduction, while increasing reduced-dimension model accuracy in cases where the zero-flux assumption does not hold. The philosophy of the method is simple, using diffusion on the membrane as an illustrative example: If it is possible to estimate the concentrations either side of the focal plane ( 𝑢 𝑖,𝐽+1𝜏 and 𝑢 𝑖,𝐽−1𝜏 ) using the concentrations on the focal plane ( 𝑢 𝑖,𝐽𝜏 ), then we can estimate the flux through the focal plane, and solve the 2D diffusion equation within the 1D reduced-dimension model. As dimension-reduced models only exist in systems exhibiting symmetry one can exploit the very symmetry used for dimension reduction to estimate concentrations on either side of the focal plane. Consider a polarity patch on the body of an elongated cell, Fig. 1A, D, F. A property of radial symmetry is that the concentration profile along all lines running through the center of a patch are identical, Fig. 1A. A property of the 1D reduced-dimension model of this system is that we know the concentration profile along a line running through the center of the patch, row 𝐽 , the focal plane. Thus, we can use interpolation on the concentration profile along row 𝐽 to estimate concentrations on rows 𝐽 ± 1 , Fig. 1F-G, Methods, SI Appendix 1. The explicit 1D finite difference diffusion equation with unrestricted movement (explicit 1D-uFDM) is, 𝑢 𝑖𝜏+1 = 𝑢 𝑖𝜏 + ∆𝑡∆𝑥 𝐷(𝑢 𝑖−1𝜏 − 2𝑢 𝑖𝜏 + 𝑢 𝑖+1𝜏 ) + ∆𝑡∆𝑦 𝐷(𝑢̃ 𝑖,𝐽−1𝜏 − 2𝑢 𝑖𝜏 + 𝑢̃ 𝑖,𝐽+1𝜏 ) where 𝑢̃ 𝑖,𝐽±1𝜏 denotes the estimated concentrations on rows 𝐽 ± 1 . For discussion on the explicit 1D-uFDM numerically stability condition see Methods. A fully implicit 1D-uFDM is ill defined, SI Appendix 2, however a semi-implicit 1D-uFDM and numerical stability condition can be derived, Methods. Illustrative example 1: Diffusion of molecules, 𝑢 , from a concentrated patch on the membrane (Fig. 2A) was modelled using the 2D diffusion equation (2D-FDM), the 1D diffusion equation (1D-FDM) and the new 1D diffusion equation with unrestricted movement (1D-uFDM), Methods. The dynamics of 𝑢 in both 1D reduced-dimension models were compared with the dynamics of 𝑢 on the focal plane of the 2D model (see SI Appendix 8 for implicit 2D-FDM, 1D-FDM and semi-implicit 1D-uFDM comparisons). Results: The reduced-dimension model using the 1D diffusion equation to describe diffusive movement through the focal plane was quantitatively inaccurate, Fig. 2B-C. On the 2D membrane molecules diffuse out of the focal plane resulting in a reduction in the mean concentration of molecules, Fig. 2D. However, molecules in the 1D-FDM model are trapped in the focal plane and the mean concentration of molecules remains constant, resulting in a higher homogeneous steady state, Fig. 2B-D. The new 1D diffusion equation with unrestricted movement (1D-uFDM) was significantly more accurate, capturing the dynamics of molecular movement through the focal plane, Fig. 2B-E. For further 1D-uFDM error analysis see SI Appendix 7. Illustrative example 2: To highlight the limitations of quantitative inaccuracy in reduced-dimension models and showcase the utility of 1D-uFDMs, we compared FRAP simulation results from the focal plane of the 2D model with the 1D reduced-dimension models. The initial fluorophore concentrations, 𝑢 , were set to reflect a uniformly covered membrane after bleaching with a Guassian laser (18), Fig. 3A. Diffusive movement was modelled using the 2D diffusion equation (2D-FDM), the 1D diffusion equation (1D-FDM) and the 1D diffusion equation with unrestricted movement (1D-uFDM). Results: For all FRAP simulations molecules diffused into the bleached area, Fig. 3B. The 1D-FDM reduced-dimension model produced a low homogeneous steady state compared to the 2D-FDM, Fig. 3B, as there was no flux of molecules into the focal plane, Fig. 3C. The 1D-FDM no flux assumption also resulted in a smaller mobile fraction than the 2D-FDM when FRAP analysis was performed on the simulated data Fig. 3D, Methods. In the 1D-FDM molecules could only enter the bleached area from a point at either side, rather than from the entire circumference of the bleached area as in the 2D model, thus the diffusion coefficient calculated from 1D-FDM data was half its actual size, Fig. 3E, Methods. The 1D-uFDM matched all 2D model 4 data, Fig. 3B-E. These results show that data obtained from membrane FRAP experiments could be used to fit parameters in 1D reduced-dimension models using 1D diffusion with unrestricted movement (1D-uFDM) but not those using the standard 1D diffusion equation. For further analysis see SI Appendix 9. Illustrative example 3: We asked to what extent the 1D reduced-dimension models could reproduce 2D reaction-diffusion (RD) dynamics. A two component, mass conserved, substrate depletion model was used to address the question (8), Fig. 4A, Methods, 𝜕 𝜕𝑡 𝑢 = 𝑢 𝑢 − 𝑢 + (0.01)∇ 𝑢 𝜕 𝜕𝑡 𝑢 = −(𝑢 𝑢 − 𝑢 ) + ∇ 𝑢 ∇ denotes the Laplacian, diffusion. Initial concentrations were set to represent a signalling event which caused a pulse conversion of molecules to 𝑢 from the 𝑢 pool, Fig. 4B, after which the dynamics of the system were described by the RD equations. The model was used to explore the effect of a progressively stronger initial signalling pulse on the dynamics of each solution. Results: Larger initial signalling pulses deplete local 𝑢 such that the RD positive feedback becomes ineffective ( 𝑢 𝑢 is very small) and a 𝑢 trough is soon formed at the location of the initial pulse, Fig. 4C, E. Recall the FRAP discussion, molecules were replenished more slowly in the reduced dimension model using standard 1D diffusion than in the 2D and 1D-uFDM models. Thus, in the 1D-FDM RD model, as the initial signalling pulse increases, creating a larger trough, movement into the trough becomes insufficient to reach the centre resulting in two narrow peaks, Fig. 4C-D. The 1D-FDM solution is qualitatively different from the 2D-FDM solution for 𝛼 ≥ 1.5 , Fig. 4C-E, SI Appendix 10. If we were unaware that the 1D-FDM solution did not represent the full 2D system then the 1D solution would be used to hypothesise a pattern on the 2D membrane. As the initial condition was a circular pulse one may hypothesise that the two peaks in the 1D solution were equivalent to a ring on the 2D membrane, giving a qualitatively inaccurate understanding of the system. The 1D diffusion equation with unrestricted movement (1D-uFDM) reproduced all the RD molecular dynamics through the focal plane, Fig. 4C-E, SI Appendix 10. Discussion
Throughout the manuscript the discussion has been focused on accurately approximating molecular movement through the focal plane. The discussion could also be viewed from a geometry-centric standpoint. It is understood that cell geometry influences molecular movement and thus cell geometry influences patterning (11, 13, 19-21) , SI Appendix 10. When one reduces the dimension of a system and restricts molecular movement to the reduced dimension they are changing the geometry of the problem. Here we show that the altered geometry of a reduced-dimension model may not approximate the geometry of the full system, leading to inaccurate results. Testing if the zero flux assumption holds would enable one to determine if the reduced-dimension geometry was a true representation of the full geometry. In spatial models dimension reduction is a commonly used, and powerful tool. While it is accepted that cell geometry affects molecular movement and hence pattern formation, reduced-dimension models are not accompanied by discussion on the relationship between molecular movement in the full and reduced dimensions. Based on the work above, we propose that such a discussion is essential to the interpretation of reduced-dimension modelling results and should be an integral part of reduced-dimension models. We have shown that the 1D diffusion equation can only be used to approximate 2D diffusion in reduced-dimension models for a specific set of systems, those for which the zero flux assumption holds. For a subset of systems that do not satisfy the zero flux assumption, we have presented a simple numerical method (the 1D-uFDM) which can be used to solve the 2D diffusion equation in a 1D reduced-dimension model by exploiting the very symmetry used to enable dimension 5 reduction. Furthermore, the methodology used to derive the 1D-uFDM is general and can be used to develop other increased accuracy reduced-dimension molecular movement equations.
Materials and Methods
Generating the interpolation mesh in 1D
Consider a 1D reduced-dimension model, reduced from a 2D uniform mesh (Fig. 1a). To solve the 1D-uFDM one must use the concentrations on the 1D mesh, row 𝐽 of the 2D mesh, to estimate concentrations at mesh-points in rows 𝐽 ± 1 of the 2D mesh (Fig. 1d, f and g). The reduced-dimension model is a ring and so has periodic boundary conditions. Without loss of generality, we will assume the centre of the patch at 𝑖 = 𝑁/2 . To generate the mesh to be interpolated from re-index the 1D mesh using 𝑛 = 𝑖 − 𝑁/2 . The transformed index, to be interpolated from is 𝑛 ={(1 − 𝑁2 ) , (2 − 𝑁2 ) , … − 1,0,1, … ( 𝑁2 − 1) , 𝑁2 } . To generate the mesh to be interpolated to use the transformed index to calculate the distance of each point in rows 𝐽 ± 1 from the centre of the patch using the equation √(n∆x) + ∆y for 𝑛 = {(1 − 𝑁2 ) , (2 − 𝑁2 ) , … − 1,0,1, … ( 𝑁2 − 1) , 𝑁2 } , Fig. 1g. Then interpolate. For systems with radial symmetry one need only estimate concentrations at a quarter of the phantom points to solve the 1D-uFDM, for example, 𝑛 = {0,1, … ( 𝑁2 − 1) , 𝑁2 } . The accuracy of estimating concentrations in rows 𝐽 ± 1 of the 2D mesh, using interpolation, is discussed in SI Appendix 1.
Explicit 1D-uFDM: derivation
The explicit 2D-FDM is, 𝑢 𝑖,𝑗𝜏+1 = 𝑢 𝑖,𝑗𝜏 + 𝑑 𝑥 (𝑢 𝑖−1,𝑗𝜏 − 2𝑢 𝑖,𝑗𝜏 + 𝑢 𝑖+1,𝑗𝜏 ) + 𝑑 𝑦 (𝑢 𝑖,𝑗−1𝜏 − 2𝑢 𝑖,𝑗𝜏 + 𝑢 𝑖,𝑗+1𝜏 ) where 𝑑 𝑥 = ∆𝑡∆𝑥 𝐷 and 𝑑 𝑦 = ∆𝑡∆𝑦 𝐷 . 𝐽 is the row of mesh points through the centre of the patch, the focal plane. When estimating 2D diffusion, through the centre of the patch, in 1D space, we only have information about row 𝐽 . Thus, we have to estimate the concentrations 𝑢 𝑖,𝐽−1𝜏 and 𝑢 𝑖,𝐽+1𝜏 . To achieve this we use the property of radial symmetry exhibited by a patch of proteins. Transform the ∆𝑥 mesh points such that the centre of the patch is at mesh point 𝑛 = 0 (see previous Methods section). The property of radial symmetry dictates that, for 𝑘 ≥ 0 , 𝑢 𝑘,𝐽+1𝜏 = 𝑢 𝑘,𝐽−1𝜏 = 𝑢 −𝑘,𝐽+1𝜏 = 𝑢 −𝑘,𝐽−1𝜏 (Fig. 1f and g). Denote the estimated concentration 𝑢 𝑛,𝐽±1𝜏 as 𝑢 √(𝑛∆𝑥) +∆𝑦 𝜏 ∀ 𝑛 = {(1 − 𝑁2 ) , … 𝑁2 } (Fig. 1g). Substituting the interpolated values of 𝑢 𝑛,𝐽±1𝜏 into the 2D-FDM equations, and removing the 𝐽 subscript we get the 1D-uFDM, 𝑢 𝑛𝜏+1 = 𝑢 𝑛𝜏 + 𝑑 𝑥 (𝑢 𝑛−1𝜏 − 2𝑢 𝑛𝜏 + 𝑢 𝑛+1𝜏 ) + 2𝑑 𝑦 (𝑢 √(𝑛∆𝑥) +∆𝑦 𝜏 − 𝑢 𝑛𝜏 ) Explicit 1D-uFDM: solution
To solve the explicit 1D-uFDM we write the 1D-uFDM equation in matrix form, 𝑢 𝜏+1 = 𝐴𝑢 𝜏 + 2𝑑 𝑦 𝑢̃ 𝐽𝜏 Where; 𝑢 𝜏 denotes the 𝑁 x vector of concentrations 𝑢 𝑛𝜏 on mesh points 𝑛 = {(1 − 𝑁2 ) , … 𝑁2 } , 𝑢̃ 𝐽𝜏 denotes the 𝑁 x vector of interpolated concentrations 𝑢 √(𝑛∆𝑥) +∆𝑦 𝜏 , 𝑛 = {(1 − 𝑁2 ) , … , 𝑁2 } , and 6 𝐴 = [ 1 − 2(𝑑 𝑥 + 𝑑 𝑦 ) 𝑑 𝑥 𝑑 𝑥 𝑥 + 𝑑 𝑦 ) 𝑑 𝑥 𝑑 𝑥 ⋱ 𝑑 𝑥 𝑑 𝑥 𝑥 + 𝑑 𝑦 ) 𝑑 𝑥 𝑑 𝑥 𝑥 + 𝑑 𝑦 )] a tridiagonal 𝑁 x 𝑁 matrix, with periodic boundary conditions. Using interpolation on the matric form we derive the solution to the explicit 1D-uFDM, 𝑢 𝜏 = 𝜆 𝜏 𝑢 + 2𝑑 𝑦 ∑ 𝜆 𝜏−𝑘 𝑢̃ 𝐽𝑘−1𝜏𝑘=1
Where 𝜆 is the vector of eigenvalues for 𝐴 . Explicit 1D-uFDM: numerical stability condition
To calculate the stability condition for 1D-uFDM numerical stability recall that the values 𝑢̃ 𝐽𝜏 are interpolated from 𝑢 𝜏 at time 𝜏 , thus, as long as a stable interpolation method is used, the solution will be stable if |𝜆| ≤ 1 ∀ 𝜆 ∈ 𝜆 . Gerschgorin’s circle theorem (17) states that 𝜆 is bounded by the inequality, 𝑥 + 𝑑 𝑦 ) − 2𝑑 𝑥 ≤ 𝜆 ≤ 1 − 2(𝑑 𝑥 + 𝑑 𝑦 ) + 2𝑑 𝑥 Which can be simplified to, 𝑥 − 2𝑑 𝑦 ≤ 𝜆 ≤ 1 − 2𝑑 𝑦 Thus, the 1D-uFDM solution will be numerically stable if, 𝑦 ≤ 1 and −1 ≤ 1 − 4𝑑 𝑥 − 2𝑑 𝑦 . 𝑦 ≤ 1 is always satisfied as 𝑑 𝑦 > 0 . The inequality, −1 ≤ 1 − 4𝑑 𝑥 − 2𝑑 𝑦 leads to the explicit 1D-uFDM stability condition, 𝑥 + 𝑑 𝑦 ≤ 1 The explicit 1D-uFDM stability condition is numerically verified in SI Appendix 3.
Semi-implicit 1D-uFDM: derivation
A fully implicit 1D-uFDM is ill defined as flux through the focal plane is inferred using the concentrations on the focal plane, SI Appendix 2. A semi-implicit numerical solver can be defined in which flux through the focal plane is solved explicitly and flux on the focal plane is solved implicitly. The semi-implicit 2D-FDM equation is, 𝑢 𝑖,𝑗𝜏+1 = 𝑢 𝑖,𝑗𝜏 + 𝑑 𝑥 (𝑢 𝑖−1,𝑗𝜏+1 − 2𝑢 𝑖,𝑗𝜏+1 + 𝑢 𝑖+1,𝑗𝜏+1 ) + 𝑑 𝑦 (𝑢 𝑖,𝑗−1𝜏 − 2𝑢 𝑖,𝑗𝜏 + 𝑢 𝑖,𝑗+1𝜏 ) Using the same reasoning used for the derivation of the explicit 1D-uFDM, the semi-implicit 1D-uFDM is, 𝑢 𝑛𝜏+1 = 𝑢 𝑛𝜏 + 𝑑 𝑥 (𝑢 𝑛−1𝜏+1 − 2𝑢 𝑛𝜏+1 + 𝑢 𝑛+1𝜏+1 ) + 2𝑑 𝑦 (𝑢 √(𝑛∆𝑥) +∆𝑦 𝜏 − 𝑢 𝑛𝜏 ) Note, the semi-implicit 1D-uFDM scheme mirrors the 1D Crank-Nicolson method which converges and is unconditionally stable(17). We will show that the same is not true for the semi-implicit 1D-uFDM. However, the semi-implicit 1D-uFDM has less strict numerical stability conditions than the explicit 1D-uFDM.
Semi-implicit 1D-uFDM: solution
To solve the semi-implicit 1D-uFDM we write it in matrix form, 𝑢 𝜏+1 = 𝐶 −1 ((1 − 2𝑑 𝑦 )𝑢 𝜏 + 2𝑑 𝑦 𝑢̃ 𝐽𝜏 )
7 where,
𝐶 = [ 1 + 2𝑑 𝑥 −𝑑 𝑥 −𝑑 𝑥 𝑥 −𝑑 𝑥 −𝑑 𝑥 ⋱ −𝑑 𝑥 −𝑑 𝑥 𝑥 −𝑑 𝑥 −𝑑 𝑥 𝑥 ] a tridiagonal 𝑁 x 𝑁 matrix, with periodic boundary conditions. Again, using iterpolation on the matric form of the equation we derive the solution to the semi-implicit 1D-uFDM, 𝑢 𝑘+1 = 1𝜆 𝑘+1 (1 − 2𝑑 𝑦 ) 𝑘+1 𝑢 + 2𝑑 𝑦 ∑(1 − 2𝑑 𝑦 ) 𝑙−1 𝑙 𝑢̃ 𝐽𝑘−𝑙𝑘𝑙=1
Where 𝜆 is the vector of eigenvalues for 𝐶 . Semi-implicit 1D-uFDM: numerical stability condition
For the solution to the semi-implicit 1D-uFDM to be stable two inequalities must hold, |1/𝜆| ≤ 1 , and |1 − 2𝑑 𝑦 | ≤ 1 . The first inequality can be investigated using Gerschgorin’s circle theorem (17). Gerschgorin’s circle theorem states that, 𝑥 − 2𝑑 𝑥 ≤ 𝜆 ≤ 1 + 2𝑑 𝑥 + 2𝑑 𝑥 Which can be simplified to, 𝑥 Thus, the inequality is always satisfied. The second inequality |1 − 2𝑑 𝑦 | ≤ 1 expands to, −1 ≤ 1 − 2𝑑 𝑦 ≤ 1 . 𝑦 ≤ 1 as 𝑑 𝑦 ≥ 0 . −1 ≤ 1 − 2𝑑 𝑦 leads to the semi-implicit 1D-uFDM stability condition, 𝑑 𝑦 ≤ 1 The semi-implicit 1D-uFDM stability condition is numerically verified in SI Appendix 3.
Illustrative example 1: Parameters
For results shown in Fig. 1i to l, to enable comparisons between the different numerical method solutions, the explicit 2D-FDM, 1D-FDM and 1D-uFDM were all solved using the same parameters:
𝐷 = 0.1 𝜇𝑚 𝑠 −1 , Δ𝑥 = Δ𝑦 = 0.1 𝜇𝑚 , Δ𝑡 = 0.01 𝑠 . These values were chosen to satisfy the numerical stability conditions of all three numerical methods. The explicit 1D-uFDM accuracy analysis found in SI Appendix 6 and 7 was also taken into consideration.
Illustrative example 2: Parameters and FRAP analysis
For the comparative FRAP analysis,
𝐷 = 0.1 𝜇𝑚 𝑠 −1 , Δ𝑥 = Δ𝑦 = 0.1 𝜇𝑚 , Δ𝑡 = 0.01 𝑠 for all solutions. The 2D ROI used to calculate the FRAP recovery curve was a circle radius and the 1D ROI a line in length (Fig. 2a). Both ROIs were placed in the centre of the bleached region. 𝑡 was calculated using linear interpolation on the FRAP recovery curve, Fig. 2d. To estimate the diffusion coefficient the equation, 𝐷̃ = 𝑟 𝑛2 + 𝑟 𝑒2 was used (18) , where 𝑟 𝑒 represents the value of the effective radius and 𝑟 𝑛 the laser radius. For the FRAP simulations, 𝑟 𝑒 = √2 and 𝑟 𝑛 = 0.55 . 8 Illustrative example 3: RD parameters and simulations
Symmetry breaking parameters, and a mass sufficient for the system to exhibit saturation due to substrate depletion ( 𝑚 = 3 ) were chosen (8). The mass of the system was defined as the mean concentration of molecules, 𝑚 = 〈𝑢 〉 + 〈𝑢 〉 where, for 𝑘 = {1,2} , 〈𝑢 𝑘 〉 = 1𝑁𝑀 ∑ ∑ 𝑢 𝑘(𝑖,𝑗)𝑀𝑗=1𝑁𝑖=1 for the 2D case, and, 〈𝑢 𝑘 〉 = 1𝑁 ∑ 𝑢 𝑘(𝑖)𝑁𝑖=1 for 1D. Δ𝑥 = Δ𝑦 = 0.1 𝜇𝑚 . Solving the RD equations was a done using a two-step process: reactions were solved using Euler’s method at time steps Δ𝑡 𝑅 = 1x10 −5 𝑠 for all simulations, diffusion was solved using explicit 2D-FDM, 1D-FDM and 1D-uFDM at time steps Δ𝑡 𝐷 = 0.002 𝑠 to ensure numerical stability of the explicit 2D-FDM (Fig. 2h to j, SI Appendix 10). Acknowledgments
This work was funded by the BBSRC, Project Grant
BB/N009339/1. Figures
Figure 1. (A)
Dimension reduction example showing a polarity patch on a spherical cell and body of an elongated cell with their reduced-dimension model representation, a stripe on a ring. White lines on 2D surfaces show the focal plane. (B)
Molecules on a 2D surface moving through the focal plane, shown as a white stripe. (C)
1D reduced-dimension model of Fig. 1D. Molecules restricted to the focal plane, shown as a white stripe. (D)
Regular mesh over a polarity patch on the body of an elongated cell, the zero flux assumption is not satisfied. (E)
Spherical mesh over a polarity patch on a spherical cell. The zero flux assumption is satisfied, 𝑢 𝑖,𝐽−1𝜏 = 𝑢 𝑖,𝐽𝜏 = 𝑢 𝑖,𝐽+1𝜏 . (F) Regular mesh over a polarity patch on the body of an elongated cell. Distances from the center of the patch of each point in rows
𝐽 ± 1 found using Pythagoras, white triangles and text. Concentrations in circles are the same because of radial symmetry. (G)
Using interpolation (black line) on the concentrations in the focal plane (orange dots) to estimate 𝑢̃ 𝑖,𝐽±1𝜏 (black open circles). J J+1
J-1 i+1 i-1 i E. B. c o n c e n t r a t i o n ( µ M ) A. C. J J+1
J-1 i+1 i-1 i D. √(2∆𝑥) + ∆𝑦 ∆𝑥 -2 -4 -5 -6 -3 -1 −√∆𝑥 + ∆𝑦 −√(3∆𝑥) + ∆𝑦 √(4∆𝑥) + ∆𝑦 G. ∆𝑦 = ∆𝑥 Δ𝑦 ∆𝑦 ∆𝑥 ∆𝑦 F. Figure 2. (A)
Initial concentrations for diffusion solutions, 𝑢 (𝑥, 0) = 𝑒 −𝑥 and inset 𝑢 (𝑥, 𝑦, 0) =𝑒 −(𝑥 +𝑦 ) . Dotted line shows the threshold for colorplots (inset and (B)). Threshold chosen to highlight differences in dynamics and homogeneous steady states in kymographs. (B) Kymographs of the concentration in the focal plane of the 2D solution, and the 1D solutions. (C)
Maximum concentration dynamics. (D)
Mean concentration in the 1D solutions compared with the mean concentration in the focal plane of the 2D solution. (E)
Mean squared distance (MSD) between the focal plane in the 2D solution and the 1D solutions. 1D-FDM comparison shown in blue, 1D-uFDM comparison shown in red. 𝜇 m initial concentrations (µM) A. time (s) m e a n c o n c e n t r a t i o n ( µ M ) D. time (s) m a x i m u m c o n c e n t r a t i o n ( µ M ) C. t i m e ( s ) c o n c e n t r a t i o n ( µ M ) 𝜇 m 𝜇 m 𝜇 m B. time (s) M S D f o r D - F D M x10 -5 M S D f o r D - u F D M E. Figure 3. (A)
Initial concentrations for FRAP simulations after bleaching with a Guassian laser, radius (18). 𝑢 (𝑥, 0) = 1 − 𝑒 −𝑥 , inset 𝑢 (𝑥, 𝑦, 0) = 1 − 𝑒 −(𝑥 +𝑦 ) . Dotted lines show the region of interest (ROI). (B) Kymographs of the concentration in the focal plane of the 2D solution, and the 1D reduced-dimension model solutions. (C)
Mean concentration in the 1D solutions compared with the mean concentration in the focal plane of the 2D solution. (D)
FRAP recovery curves. The 2D ROI is a circle diameter , 1D a line length . (E) Half time, 𝑡 , actual diffusion coefficients, 𝐷 , and diffusion coefficients estimated using simulated data, 𝐷̃ . 12 Figure 4. (A)
Cartoon of the RD model. Black arrows represent reactions, green dashed arrows represent diffusion. Mass is conserved, molecules are only changed, never created or lost. Substrate 𝑢 is depleted in the creation of 𝑢 . (B) Initial concentrations of 𝑢 and 𝑢 for RD simulations, 𝑢 (𝑥, 𝑦, 0) = 𝛼𝑒 −(𝑥 +𝑦 ) , 𝑢 (𝑥, 0) = 𝛼𝑒 −𝑥 , and 𝑢 = 𝑚 − 𝑢 , 𝛼 = [0.5, 1, … , 3] . (C) Kymographs of the 𝑢 concentration in the focal plane of the 2D solution and the 1D solutions. (D) Steady state colorplots of the 2D and 1D 𝑢 solutions. (E) Central 𝑢 concentration in the 1D and 2D solutions. Describing the movement of molecules in reduced-dimension models.
Natasha S. Savage. University of Liverpool, Liverpool, UK.
Email: [email protected]
Supplementary Information Appendix 1: Accuracy of estimating concentrations at phantom points using interpolation.
Numerical solutions within this Brief Communication were performed in MATLAB. Spline interpolation was found to be the most accurate interpolation method in MATLAB for the concentration profiles investigated here (data not shown). However, if a concentration profile approaches zero with a steep gradient spline interpolation can return a small negative concentration. Any negative concentrations resulting from MATLAB spline interpolation were set to zero. To investigate the accuracy of estimating concentrations at phantom mesh-points (those in rows
𝐽 ± 1 ) using MATLAB spline interpolation, 2D concentration profiles were generated and the actual concentrations on rows
𝐽 ± 1 were compared with the estimated concentrations, which had been interpolated from the 1D concentration profile (Methods). To quantify interpolation accuracy the sum squared distance (SSD) between the interpolated and actual concentration profiles was calculated. The concentration profiles used to investigate the accuracy of estimating concentrations at phantom points can be seen in Fig. S1. To better understand the influence of ∆𝑥 and ∆𝑦 on the estimation of phantom Point concentrations, a range of ∆𝑥 and ∆𝑦 values, between and , was used for each concentration profile. As expected, increasing ∆𝑥 decreases the accuracy of the interpolation (Fig. S2, column 1). Comparing columns 1 and 2 of Fig. S2 we see that the error around the concentration peak, in the central 10% of the mesh, accounts for most of the interpolation error. Indeed, the greater the curvature of the concentration profile the greater the interpolation error. In order to estimate concentrations at phantom points on the special boundaries one must interpolate outside the data (Fig. S5b). Fig. S2, column 3, shows that, as we may expect, the greater ∆𝑦 the greater the error in estimation outside the boundary. Supplementary Table 1 shows the maximum and minimum SSD when estimating the phantom points of all the concentration profiles shown in Fig. S1a, c, and e. The estimation error is small for small ∆𝑥 and ∆𝑦 values (Fig S2 and Table S1). Decreasing ∆𝑥 increases computation time, to reduce estimation error for larger ∆𝑥 a small ∆𝑦 should be chosen (Fig. S2). The choices of ∆𝑥 and ∆𝑦 should be chosen to satisfy the stability conditions of the numerical method (Methods). Appendix 2: An implicit 1D-uFDM is ill defined.
The implicit 2D-FDM is, 𝑢 𝑖,𝑗𝜏+1 = 𝑢 𝑖,𝑗𝜏 + ∆𝑡∆𝑥 𝐷(𝑢 𝑖−1,𝑗𝜏+1 − 2𝑢 𝑖,𝑗𝜏+1 + 𝑢 𝑖+1,𝑗𝜏+1 ) + ∆𝑡∆𝑦 𝐷(𝑢 𝑖,𝑗−1𝜏+1 − 2𝑢 𝑖,𝑗𝜏+1 + 𝑢 𝑖,𝑗+1𝜏+1 ) Using the same reasoning as for the derivation of the explicit 1D-uFDM, the implicit 1D-uFDM is, 14 𝑢 𝑛𝜏+1 = 𝑢 𝑛𝜏 + ∆𝑡∆𝑥 𝐷(𝑢 𝑛−1𝜏+1 − 2𝑢 𝑛𝜏+1 + 𝑢 𝑛+1𝜏+1 ) + 2 ∆𝑡∆𝑦 𝐷 (𝑢 √(𝑛∆𝑥) +∆𝑦 𝜏+1 − 𝑢 𝑛𝜏+1 ) The vector form of the implicit 1D-uFDM is, 𝑢 𝜏 = 𝐵𝑢 𝜏+1 + 2𝑑 𝑦 𝑢̃ 𝐽𝜏+1
𝐵 = [ 1 + 2(𝑑 𝑥 + 𝑑 𝑦 ) −𝑑 𝑥 −𝑑 𝑥 𝑥 + 𝑑 𝑦 ) −𝑑 𝑥 ⋱ −𝑑 𝑥 𝑥 + 𝑑 𝑦 ) −𝑑 𝑥 −𝑑 𝑥 𝑥 + 𝑑 𝑦 )] a tridiagonal 𝑁 x 𝑁 matrix, with periodic boundary conditions. The solution to the implicit 1D-uFDM is, 𝑢 𝜏+1 = 𝐵 −1 (𝑢 𝜏 − 2𝑑 𝑦 𝑢̃ 𝐽𝜏+1 ) However, the values of 𝑢 𝜏+1 needed to interpolate the values 𝑢̃ 𝐽𝜏+1 , and so 𝑢̃ 𝐽𝜏+1 in the above equation are not defined.
Appendix 3: Numerical verification of stability Conditions.
The explicit and semi-implicit 1D-uFDM numerical stability conditions were tested numerically. Recall, 𝐷 ∆𝑡∆𝑥 = 𝑑 𝑥 and 𝐷 ∆𝑡∆𝑦 = 𝑑 𝑦 . Consider first the explicit 1D-uFDM stability condition in terms of ∆𝑥 , ∆𝑦 , ∆𝑡 and 𝐷 , ∆𝑡∆𝑥 + 𝐷 ∆𝑡∆𝑦 ≤ 1 . For 𝛽 ≥ 0 we would expect a numerical solution with ∆𝑡∆𝑥 +𝐷 ∆𝑡∆𝑦 = 1 + 𝛽 to be numerically unstable and a numerical solution with ∆𝑡∆𝑥 + 𝐷 ∆𝑡∆𝑦 = 1 − 𝛽 to be numerically stable. To test the explicit 1D-uFDM stability condition ∆𝑥 , ∆𝑦 and 𝐷 were set and ∆𝑡 ± was calculated using the formula,
2𝐷 ∆𝑡 ± ∆𝑥 + 𝐷 ∆𝑡 ± ∆𝑦 = 1 ± 𝛽 ∆𝑡 ± = 1 ± 𝛽D ( ∆𝑥 ∆𝑦 + ∆𝑥 ) Numerical solutions were run for ∆𝑥 = 0.01,0.02, … ,0.1 𝜇𝑚 , ∆𝑦 = 0.01,0.02, … ,0.1 𝜇𝑚 , 𝐷 =0.1 𝜇𝑚 𝑠 −1 and ∆𝑡 = ∆𝑡 ± 𝑠 , 𝛽 = 1x10 −4 , 1x10 −3 , 1x10 −2 . The initial condition for each simulation was 𝑢(𝑥, 𝑡) = 𝑢(𝑥, 0) = 𝑒 −𝑥 for 𝑥 ∈ [0,12] 𝜇𝑚 (Fig. S1a, 𝛼 = 1 ). Numerical solutions were calculated for one hour, or until they crashed as a result of numerical instability. Fig. S3a shows the stability boundaries being tested. Table S2 shows the smallest values of 𝛽 , out of those tested, for which numerical solutions with ∆𝑡 + were numerically unstable and ∆𝑡 − were numerically stable. For the values of ∆𝑥, ∆𝑦 in the blue cells of Table S2, ∆𝑡 − simulations with 𝛽 = 1x10 −4 , 1x10 −3 were numerically unstable. To understand how far from the stability boundary the blue celled values in Table S2 were, the ∆𝑡 ± equation was rearranged such that ∆𝑡 ± = Δ𝑡 ± 𝛽Δ𝑡 , where ∆𝑡 = ( ∆𝑥 ∆𝑦 +∆𝑥 ) . The values Δ𝑡 ± (1x10 −2 )Δ𝑡 are shown in Table S3. The semi-implicit 1D-uFDM stability condition was tested numerically using the same methodology, with ∆𝑡 ± = ∆𝑦 (Fig. S3b, Table S4). Simulations with ∆𝑡 + were numerically unstable and ∆𝑡 − were numerically stable for all ∆𝑥 , ∆𝑦 simulated and 𝛽 = 1x10 −4 . 15 Appendix 4: 1D Fourier solution.
Regarding notation: For the finite difference scheme, 𝑢 𝑛𝜏 denotes the concentration of a species on mesh point 𝑛 at time 𝜏 . To avoid confusion, here we define the notation for the Fourier solution. Let 𝑣(𝑥, 𝑡) or 𝑣(𝑥, 𝑦, 𝑡) describe the concentration of a species at position 𝑥 or (𝑥, 𝑦) at time 𝑡 . The Fourier general solution will be comprised of the sum of fundamental solutions, using the principle of superposition. 𝑣 𝑛 (𝑥, 𝑡) or 𝑣 𝑛 (𝑥, 𝑦, 𝑡) will denote the 𝑛 th fundamental solution. For this investigation, the 1D initial condition was 𝑣(𝑥, 0) = 𝑒 −𝑥 , 2D 𝑣(𝑥, 𝑦, 0) = 𝑒 −(𝑥 +𝑦 ) (Fig. S1a and b, 𝛼 = 1 ), and 𝐷 = 0.1 𝜇𝑚 𝑠 −1 . The 1D Fourier solution informs the 2D solution so we will present that first. Consider a 1D space – 𝐿 to 𝐿 , with periodic boundary conditions. Diffusive movement of 𝑣(𝑥, 𝑡) in the 1D space is given by the equation, 𝜕 𝜕𝑡 𝑣(𝑥, 𝑡) = 𝐷 𝜕 𝜕𝑥 𝑣(𝑥, 𝑡) , where 𝐷 is the diffusion coefficient. The periodic boundary conditions are, 𝑣(−𝐿, 𝑡) = 𝑣(𝐿, 𝑡) Boundary condition 1 𝜕 𝜕𝑥 𝑣(𝑥, 𝑡)| −𝐿 = 𝜕 𝜕𝑥 𝑣(𝑥, 𝑡)| 𝐿 Boundary condition 2
The diffusion equation can be solved by separation of variables and Fourier series expansions (22, 23). Separation of variables states 𝑣(𝑥, 𝑡) will have the form, 𝑣(𝑥, 𝑡) = 𝑋(𝑥)𝑇(𝑡) . Substitute 𝑣(𝑥, 𝑡) =𝑋(𝑥)𝑇(𝑡) into the diffusion equation to get,
1𝐷 𝑇(𝑡) 𝜕 𝜕𝑡
𝑇(𝑡) = 𝜕𝑥 𝑋(𝑥)
Equation 3
Both sides of the Equation 3 must equal a constant as one side depends only on time and the other depends only on space. Let the constant be −𝜆 . The time aspect of Equation 3 is, 𝜕 𝜕𝑡 𝑇(𝑡) +𝜆 𝐷 𝑇(𝑡) = 0 , which has the solution,
𝑇(𝑡) = 𝐶𝑒 −𝜆𝐷𝑡 , where 𝐶 is an arbitrary constant. The space aspect of Equation 3 is, 𝜕 𝜕𝑥 𝑋(𝑥) + 𝜆 𝑋(𝑥) = 0
Equation 4
There are three solutions to Equation 4, one for 𝜆 < 0, 𝜆 = 0 and 𝜆 > 0 . The validity of these solutions are tested using boundary conditions 1 and 2. For 𝜆 < 0 , 𝑋(𝑥) = 𝑘 cosh(√−𝜆 𝑥) + 𝑘 sinh(√−𝜆 𝑥) , which has only trivial solutions. For 𝜆 = 0 , 𝑋(𝑥) = 𝑘 𝑥 + 𝑘 . Applying Boundary Conditions 1 and 2 we get the fundamental solution, 𝑣 (𝑥, 𝑡) = 𝐶 Equation 5
16 For 𝜆 > 0 , Equation 4 has the solution,
𝑋(𝑥) = 𝑘 cos(√𝜆 𝑥) + 𝑘 sin(√𝜆 𝑥) . Applying Boundary Condition 1 gives, 𝑘 cos(−√𝜆 𝐿) + 𝑘 sin(−√𝜆 𝐿) = 𝑘 cos(√𝜆 𝐿) + 𝑘 sin(√𝜆 𝐿) . As cos is even we get, − 𝑘 sin(√𝜆 𝐿) = 𝑘 sin(√𝜆 𝐿) . This equality holds if 𝑘 = 0 , or 𝜆 𝑛 = ( 𝑛𝜋𝐿 ) , 𝑛 = 1,2,3 … . The solution 𝑋(𝑥) = 𝑘 cos(√𝜆 𝑥) + 𝑘 sin(√𝜆 𝑥) becomes 𝑋 𝑛 (𝑥) = 𝑘 cos ( 𝑛𝜋𝐿 𝑥) + 𝑘 sin ( 𝑛𝜋𝐿 𝑥) , which satisfies Boundary Condition 2. Note the solution for 𝑘 = 0 is accounted for in 𝑋 𝑛 (𝑥) because of the principle of superposition. Thus we have the set of fundamental solutions, 𝑣 𝑛 (𝑥, 𝑡) = (𝐴 𝑛 cos ( 𝑛𝜋𝐿 𝑥) + 𝐵 𝑛 sin ( 𝑛𝜋𝐿 𝑥)) 𝑒 −( 𝑛𝜋𝐿 ) 𝐷𝑡 Equation 6 where 𝐴 𝑛 = 𝑘 𝐶 𝑛 and 𝐵 𝑛 = 𝑘 𝐶 𝑛 . Using the principle of superposition, we write the general solution to the diffusion equation with Boundary Conditions 1 and 2 as a linear composition of all solutions (Equation 5 and Equation 6), 𝑣(𝑥, 𝑡) = 𝑣 (𝑥, 𝑡) + ∑ 𝑣 𝑛 (𝑥, 𝑡) 𝑛=1,2… , 𝑣(𝑥, 𝑡) = 𝐶 + ∑ [(𝐴 𝑛 cos ( 𝑛𝜋𝐿 𝑥) + 𝐵 𝑛 sin ( 𝑛𝜋𝐿 𝑥)) 𝑒 −( 𝑛𝜋𝐿 ) 𝐷𝑡 ] 𝑛=1,2… Equation 7
We use the initial condition and Fourier analysis to determine the coefficients of Equation 7, and thus obtain the solution of the diffusion equation with the given initial condition. The initial condition is (𝑥, 0) = 𝑒 −𝑥 . At 𝑡 = 0 Equation 7 becomes a Fourier series, 𝑒 −𝑥 = 𝐶 + ∑ [𝐴 𝑛 cos ( 𝑛𝜋𝐿 𝑥) + 𝑛=1,2… 𝐵 𝑛 sin ( 𝑛𝜋𝐿 𝑥)] . The initial condition is an even function and so 𝐵 𝑛 = 0 ∀ 𝑛 . Fourier analysis states that 𝐶 = ∫ 𝑒 −𝑥 𝐿−𝐿 𝑑𝑥 , giving, and 𝐴 𝑛 = ∫ cos ( 𝑛𝜋𝐿 𝑥) 𝑒 −𝑥 𝐿−𝐿 𝑑𝑥 = . Thus, the general solution for the 1D diffusion equation with initial condition 𝑣(𝑥, 0) = 𝑒 −𝑥 and periodic boundary conditions is, 𝑣(𝑥, 𝑡) = 𝐶 + ∑ [𝐴 𝑛 cos ( 𝑛𝜋𝐿 𝑥) 𝑒 −( 𝑛𝜋𝐿 ) 𝐷𝑡 ] 𝑛=1,2… 𝐶 = √𝜋 erf (𝐿) 𝐴 𝑛 = 𝑖√𝜋 𝑒 −( 𝑛𝜋2𝐿 ) (erfi ( 𝑛𝜋2𝐿 − 𝑖𝐿) − erfi ( 𝑛𝜋2𝐿 + 𝑖𝐿)) Appendix 5: 2D Fourier solution.
Consider a 2D space 𝑥 ∈ [−𝐿, 𝐿] and 𝑦 ∈ [−𝐿, 𝐿] with periodic boundary conditions. The diffusion equation in 2D space is 𝜕 𝜕𝑡 𝑣(𝑥, 𝑦, 𝑡) = 𝐷 ( 𝜕 𝜕𝑥 𝑣(𝑥, 𝑦, 𝑡) + 𝜕 𝜕𝑦 𝑣(𝑥, 𝑦, 𝑡)) , with periodic boundary conditions, 17 𝑣(−𝐿, 𝑦, 𝑡) = 𝑣(𝐿, 𝑦, 𝑡) 𝑣(𝑥, −𝐿, 𝑡) = 𝑣(𝑥, 𝐿, 𝑡) 𝜕 𝜕𝑥 𝑣(𝑥, 𝑦, 𝑡)| −𝐿 = 𝜕 𝜕𝑥 𝑣(𝑥, 𝑦, 𝑡)| 𝐿 𝜕 𝜕𝑦 𝑣(𝑥, 𝑦, 𝑡)| −𝐿 = 𝜕 𝜕𝑦 𝑣(𝑥, 𝑦, 𝑡)| 𝐿 Separation of variables states 𝑣(𝑥, 𝑦, 𝑡) will have the form, 𝑣(𝑥, 𝑦, 𝑡) = 𝑋(𝑥)𝑌(𝑦)𝑇(𝑡) . Substitute 𝑣(𝑥, 𝑦, 𝑡) = 𝑋(𝑥)𝑌(𝑦)𝑇(𝑡) into the 2D diffusion equation to get,
1𝐷 𝑇(𝑡) 𝜕 𝜕𝑡
𝑇(𝑡) = 𝜕𝑥 𝑋(𝑥) + 𝜕𝑦 𝑌(𝑦) = −𝜆 = −(𝜆 𝑛 + 𝜆 𝑚 ) . 𝑇(𝑡) and
𝑋(𝑥) have the same solution as they did when considering the 1D diffusion equation, and the solution for
𝑌(𝑦) will have the same form as
𝑋(𝑥) . However, for the 2D case 𝜆 𝑛 = ( 𝑛𝜋𝐿 ) , 𝜆 𝑚 = ( 𝑚𝜋𝐿 ) and 𝜆 = ( 𝑛𝜋𝐿 ) + ( 𝑚𝜋𝐿 ) . 𝑇(𝑡) = 𝐶𝑒 −(( 𝑛𝜋𝐿 ) +( 𝑚𝜋𝐿 ) )𝐷𝑡 𝑋(𝑥) = 𝑘 𝑋 𝑛 (𝑥) = 𝑘 cos ( 𝑛𝜋𝐿 𝑥) + 𝑘 sin ( 𝑛𝜋𝐿 𝑥) 𝑌(𝑦) = 𝑘 𝑌 𝑛 (𝑦) = 𝑘 cos ( 𝑚𝜋𝐿 𝑦) + 𝑘 sin ( 𝑚𝜋𝐿 𝑦) The corresponding solutions for 𝑣(𝑥, 𝑦, 𝑡) are 𝑣 𝑛𝑚 (𝑥, 𝑦, 𝑡) , 𝑣 (𝑥, 𝑦, 𝑡) = 𝐶 𝑣 (𝑥, 𝑦, 𝑡) = 𝐶 (𝐴 𝑚 cos ( 𝑚𝜋𝐿 𝑦) + 𝐵 𝑚 sin ( 𝑚𝜋𝐿 𝑦)) 𝑒 −( 𝑚𝜋𝐿 ) 𝐷𝑡 𝑣 𝑛0 (𝑥, 𝑦, 𝑡) = 𝐶 𝑛0 (𝐴 𝑛 cos ( 𝑛𝜋𝐿 𝑥) + 𝐵 𝑛 sin ( 𝑛𝜋𝐿 𝑥)) 𝑒 −( 𝑛𝜋𝐿 ) 𝐷𝑡 𝑣 𝑛𝑚 (𝑥, 𝑦, 𝑡) = (𝐴 𝑛 cos ( 𝑛𝜋𝐿 𝑥) + 𝐵 𝑛 sin ( 𝑛𝜋𝐿 𝑥)) (𝐴 𝑚 cos ( 𝑚𝜋𝐿 𝑦) + 𝐵 𝑚 sin ( 𝑚𝜋𝐿 𝑦)) 𝑒 −(( 𝑛𝜋𝐿 ) +( 𝑚𝜋𝐿 ) )𝐷𝑡
18 Using the principle of superposition, we write the general solution to the 2D diffusion equation with periodic boundary conditions, 𝑣(𝑥, 𝑦, 𝑡) = 𝑣 (𝑥, 𝑦, 𝑡) + ∑ 𝑣 (𝑥, 𝑦, 𝑡) 𝑚=1,2… + ∑ 𝑣 𝑛0 (𝑥, 𝑦, 𝑡) 𝑛=1,2… + ∑ 𝑣 𝑛𝑚 (𝑥, 𝑦, 𝑡) 𝑛,𝑚=1,2… Again, we use Fourier analysis find the Fourier coefficients and solve the general solution for a given initial condition. The general solution for the 2D diffusion equation with initial condition 𝑣(𝑥, 𝑦, 0) = 𝑒 −(𝑥 +𝑦 ) is, 𝑣(𝑥, 𝑦, 𝑡) = 𝐶 + ∑ 𝐴 cos ( 𝑚𝜋𝐿 𝑦) 𝑚=1,2… 𝑒 −( 𝑚𝜋𝐿 ) 𝐷𝑡 + ∑ 𝐴 𝑛0 cos ( 𝑛𝜋𝐿 𝑥) 𝑒 −( 𝑛𝜋𝐿 ) 𝐷𝑡𝑛=1,2… + ∑ 𝐴 𝑛𝑚 cos ( 𝑛𝜋𝐿 𝑥) cos ( 𝑚𝜋𝐿 𝑦) 𝑛,𝑚=1,2… 𝑒 −(( 𝑛𝜋𝐿 ) +( 𝑚𝜋𝐿 ) )𝐷𝑡 𝐶 = 𝜋(2𝐿) erf(𝐿) 𝐴 = −𝜋𝑖(2𝐿) erf(𝐿) 𝑒 −( 𝑚𝜋2𝐿 ) (erfi ( 𝐿𝑖 − 𝑚𝜋2𝐿 ) + erfi ( 𝐿𝑖 + 𝑚𝜋2𝐿 )) 𝐴 𝑛0 = −𝜋𝑖(2𝐿) erf(𝐿) 𝑒 −( 𝑛𝜋2𝐿 ) (erfi ( 𝐿𝑖 − 𝑛𝜋2𝐿 ) + erfi ( 𝐿𝑖 + 𝑛𝜋2𝐿 )) 𝐴 𝑛𝑚 = −𝜋(2𝐿) 𝑒 −( 𝑛𝜋2𝐿 ) −( 𝑚𝜋2𝐿 ) (erfi ( 𝐿𝑖 − 𝑛𝜋2𝐿 ) + erfi ( 𝐿𝑖 + 𝑛𝜋2𝐿 )) (erfi ( 𝐿𝑖 − 𝑚𝜋2𝐿 ) + erfi ( 𝐿𝑖 + 𝑚𝜋2𝐿 )) Appendix 6: Steady state accuracy of 1D-uFDMs.
The accuracy of the 1D-uFDM solutions was investigated numerically by comparing it with the central slice of the 2D Fourier solution (Appendix 5). The comparisons were made using three sets of in-silico experiments. In set 1, we asked how 1D-uFDM accuracy was affected by mesh size by fixing ∆𝑡 and varying ∆𝑥, ∆𝑦 . In set 2 we asked how accuracy was affected by increasing ∆𝑡 . Here, we fixed ∆𝑥, ∆𝑦 and varied ∆𝑡 . Finally, in set 3, we looked at the interdependence of ∆𝑥, ∆𝑦 and ∆𝑡 by varying all three. First we looked at the accuracy of the 1D-uFDMs when predicting steady state on the central slice of the 2D Fourier solution. Simulation set 1: To ascertain how the accuracy of the 1D-uFDMs was affected by mesh size, for each ∆𝑥 = 0.01,0.02, … ,0.1 𝜇𝑚 , a simulation was run for ∆𝑦 =0.01,0.02, … ,0.1 𝜇𝑚 . To ensure numerical stability for all ∆𝑥, ∆𝑦 , and for both explicit and semi-implicit numerical methods, ∆𝑡 was set using the explicit 1D-uFDM stability condition (Fig. S3), ∆𝑡 − (Δ𝑥, Δ𝑦, 𝐷, 𝛽) = ∆𝑡 − (0.01,0.01,0.1, 0.1) = ( ∆𝑥 ∆𝑦 +∆𝑥 ) = 0.0003 𝑠 . The initial conditions were, 𝑢 (𝑥, 𝑡) = 𝑢 (𝑥, 0) = 𝑒 −𝑥 , 𝑢 (𝑥, 𝑦, 𝑡) = 𝑢 (𝑥, 𝑦, 0) = 𝑒 −(𝑥 +𝑦 ) (Fig. S1a and b, 𝛼 = 1 ). The protocol was thus, for each ∆𝑥 , ∆𝑦 combination: The explicit and semi-implicit 1D-uFDMs were run to
300 𝑠 . The 2D Fourier solution was calculated at and
300 𝑠 . The mean squared distance (MSD) was calculated between each point on the 1D mesh and the corresponding point on the central row of the 2D mesh, for both and
300 𝑠 . If the MSD was zero at then the MSD at
300 𝑠 (steady state) was recorded. Changing mesh size did not dramatically affect the accuracy of the unrestricted finite difference methods, Fig. S4a. However, increasing ∆𝑥 resulted in a small 19 increase in accuracy and increasing ∆𝑦 resulted in a small decrease. For ∆𝑡 − = 0.0003 𝑠 the semi-implicit unrestricted movement finite difference method was marginally more accurate than the explicit method. Simulation set 2: Next, we investigated the effect of ∆𝑡 on the accuracy of the unrestricted movement finite difference methods by fixing ∆𝑥 = Δ𝑦 = 0.1 𝜇𝑚 and varying ∆𝑡 . ∆𝑥 = Δ𝑦 = 0.1 𝜇𝑚 was chosen to give the largest numerically stable range of ∆𝑡 values . For both semi-implicit and explicit 1D-uFDMs ∆𝑡 = 0.0003 to at intervals of 0.0027 s. As the semi-implicit 1D-uFDM is numerically stable for greater values of ∆𝑡 (Fig. S3b) the semi-implicit 1D-uFDM was also run for ∆𝑡 = 0.04, 0.05, … 0.1 s (Fig. S4b). For the explicit 1D-uFDM increasing ∆𝑡 had no significant effect on the accuracy. The semi-implicit 1D-uFDM’s accuracy was significantly increased with increasing ∆𝑡 . Given the difference in accuracy between the accuracy of semi-implicit and explicit 1D-FDMs for increasing ∆𝑡 one would expect that fixing ∆𝑡 > ∆𝑥 ≥ 0.05 𝜇𝑚, Δ𝑦 ≥ 0.05 𝜇𝑚 and ∆𝑡 − (0.5,0.5,0.1, 0.1) = ( ∆𝑥 ∆𝑦 +∆𝑥 ) = 0.0075 𝑠 . The results of these simulations confirmed expectations (Fig. S4c). Fixing ∆𝑥 = ∆𝑦 = 0.05 𝜇𝑚 and varying ∆𝑡 gave further confirmation of the accuracy trends reported in Fig. S4d. For the explicit 1D-uFDM the choice of ∆𝑥, ∆𝑦 and ∆𝑡 has little effect on the accuracy of the numerical method, as long as the choice satisfices the explicit 1D-uFDM stability condition. For the semi-implicit 1D-uFDM increasing ∆𝑡 benefits accuracy, however to enable larger ∆𝑡 one must increase ∆𝑦 to ensure numerical stability. Simulation set 3 explores the interdependence of ∆𝑥, ∆𝑦 and ∆𝑡 with regards to steady state accuracy. For these simulations, all three parameters were changed. ∆𝑥 = 0.01,0.02, … ,0.1 𝜇𝑚, ∆𝑦 = 0.01,0.02, … ,0.1 𝜇𝑚 , and ∆𝑡 was set for each ∆𝑥 , ∆𝑦 combination using the equation ∆𝑡 − (Δ𝑥, Δ𝑦, 𝐷, 𝛽) = ∆𝑡 − (Δ𝑥, Δ𝑦, 0.1, 0.1) = ( ∆𝑥 ∆𝑦 +∆𝑥 ) 𝑠 for explicit 1D-uFDM simulations and ∆𝑡 − (Δ𝑦, 𝐷, 𝛽) = ∆𝑡 − (Δ𝑦, 0.1, 0.1) = ∆𝑦 for semi-implicit. As expected, the explicit 1D-uFDM had the same accuracy with the larger values of ∆𝑡 as it did with ∆𝑡 = 0.0003 s (compare solid lines on Fig. S4a and S4e). As the stability boundary of the semi-implicit 1D-uFDM depends on ∆𝑦 only, increasing ∆𝑦 increases the value ∆𝑡 − (Δ𝑦, 0.1, 0.1) and thus increases the accuracy of the semi-implicit numerical method (Fig. S4e). Continually increasing the value of ∆𝑦 , and thus ∆𝑡 , should eventually reduce the accuracy of the numerical method. To investigate this ∆𝑥 was set to , ∆𝑦 was increased from to , and ∆𝑡 − (Δ𝑦, 0.1, 0.1) = ∆𝑦 𝑠 . Indeed, the minimum MSD was found to be −9 for ∆𝑦 = 0.38 𝜇𝑚 , ∆𝑡 = 1.2996 𝑠 (Fig. S4f). Taken together these investigations show that to improve the accuracy of the steady state estimation when using the explicit 1D-uFDM one should choose a larger ∆𝑥 and smaller ∆𝑦 , then set ∆𝑡 with the formula ∆𝑡 − (Δ𝑥, Δ𝑦, 𝐷, 𝛽) = ( ∆𝑥 ∆𝑦 +∆𝑥 ) . For the semi-implicit method one should choose a larger ∆𝑦 and set ∆𝑡 using the formula ∆𝑡 − (Δ𝑦, 𝐷, 𝛽) = ∆𝑦 . Appendix 7: Accuracy dynamics of 1D-uFDMs.
Next, we compare the dynamics of the 1D-uFDMs to see how faithfully they estimate the approach to steady state. Similar to the steady state accuracy investigation we performed three sets of simulations. Simulation set 1: Investigate the time course dependence of 1D-uFDMs accuracy, and its relationship to mesh size. As calculating the Fourier solution every second for
300 𝑠 is computationally time consuming we chose to investigate fewer mesh sizes than in the steady state accuracy investigation. For each ∆𝑥 = 0.02,0.04, … ,0.1 𝜇𝑚 , a simulation was run for ∆𝑦 = 𝜇𝑚 . Again, to ensure numerical stability for all ∆𝑥, ∆𝑦 , and for both explicit and semi-implicit numerical methods, ∆𝑡 was set using the explicit 1D-uFDM stability condition (Fig. S3). ∆𝑡 − (Δ𝑥, Δ𝑦, 𝐷, 𝛽) = ∆𝑡 − (0.01,0.01,0.1, 0.1) = 0.0003 𝑠 . The initial conditions were 𝑢 (𝑥, 𝑡) =𝑢 (𝑥, 0) = 𝑒 −𝑥 , 𝑢(𝑥, 𝑦, 𝑡) = 𝑢(𝑥, 𝑦, 0) = 𝑒 −(𝑥 +𝑦 ) (Fig. S1a and b, 𝛼 = 1 ). The protocol was thus, for each ∆𝑥 , ∆𝑦 combination: The explicit and semi-implicit 1D-uFDMs were run to
300 𝑠 . As ∆𝑡 did not divide exactly into one, the 1D-uFDMs concentration profiles were recorded at close to one second intervals, and the times of data collection were recorded. The 2D Fourier solutions were calculated at the same times as the 1D-uFDMs’ data was collected. For every time point the MSD was calculated between the 1D-uFDMs and the central row of the 2D Fourier solution. Fig. S5a shows the MSD time evolution. As with the steady state MSD analysis, MSD dynamics were similar for all ∆𝑥 , ∆𝑦 tested, for both the explicit and semi-implicit 1D-uFDMs. In all cases the MSD increased rapidly at around
30 𝑠 . To understand the cause of the rapid MSD increase we calculated the absolute distance between each point of the semi-implicit 1D-uFDM and the central row of the 2D Fourier solution, at every second, for
Δ𝑥 = 0.1 𝜇𝑚 , Δ𝑦 = 0.02 𝜇𝑚 as this had the smallest MSD at steady state (Fig. S4a, Fig. S5c). The absolute distance analysis showed that the absolute distance began to increase at the boundary edges and propagated through the mesh. Looking at the semi-implicit 1D-uFDM concentration profile we see that the increase in absolute distance at the boundaries corresponds to the rate of change of concentration at the boundaries (compare Fig. S5c and S5d, see also S5e). This is unsurprising as estimating concentrations 𝑢 √(𝑁/2 ∆𝑥) +∆𝑦 𝜏 , which would lie at the boundaries of the 2D domain (black circle Fig. S5b), involves a 1D interpolation into space without any concentration information (orange circle Fig. S5b). As the curvature at the boundary increases so too does the interpolation error. Looking more closely at the MSD time course data (Fig. S5a), we see a small peak before the edge error is detected (Fig. S6a and S6b). This initial MSD peak has a strong ∆𝑦 dependency. To ask if this initial MSD peak was a result of interpolation error we looked at the interpolation accuracy data (Fig. S2, row α = 1 , ∆𝑥 = 0.1 𝜇𝑚 ). This data showed that the interpolation accuracy does decrease with increasing ∆𝑦 , and that the decrease is predominantly generated in the central region, around the location of the peak. However the size of the initial MSD peak is not accounted for by the interpolation inaccuracy on the initial concentration (Fig. S5a, row 2, column 1, α = 1 , ∆𝑥 = 0.1 𝜇𝑚 , L = 12 𝜇𝑚 , divide max SSD in figure by
𝐿/∆𝑥 = 120 to get a maximum MSD around 10 -13 ). The 1D-uFDMs are composed of two parts, the 1D-FDM (highlighted in blue text in the equations below), and the terms utilising the interpolated concentrations. 𝑢 𝑛𝜏+1 = 𝑢 𝑛𝜏 + 𝑑 𝑥 (𝑢 𝑛−1𝜏 − 2𝑢 𝑛𝜏 + 𝑢 𝑛+1𝜏 ) + 2𝑑 𝑦 (𝑢 √(𝑛∆𝑥) +∆𝑦 𝜏 − 𝑢 𝑛𝜏 ) 𝑢 𝑛𝜏+1 = 𝑢 𝑛𝜏 + 𝑑 𝑥 (𝑢 𝑛−1𝜏+1 − 2𝑢 𝑛𝜏+1 + 𝑢 𝑛+1𝜏+1 ) + 2𝑑 𝑦 (𝑢 √(𝑛∆𝑥) +∆𝑦 𝜏 − 𝑢 𝑛𝜏 ) To ask if the error inherent in the 1D-FDMs could account for the initial MSD peak the 1D-uFDMs, we compared the solutions of the 1D-FDMs with the 1D Fourier solutions for the diffusion equation (1D Fourier solution derived in Appendix 4). To enable direct error comparisons between the 1D-FDMs inherent error analysis and the 1D-uFDMs error analysis the 1D-FDMs were run for
300 𝑠 , with ∆𝑥 = 0.02,0.04, … ,0.1 𝜇𝑚 , ∆𝑡 = 0.0003 𝑠 and initial condition 𝑢 (𝑥, 𝑡) = 𝑢 (𝑥, 0) = 𝑒 −𝑥 . Again, as ∆𝑡 does not divide exactly into one, the 1D-FDMs concentration profiles were recorded at close to one second intervals, and the times of data collection were recorded. The 1D Fourier solutions were calculated at the same times. For every time point, the MSDs and SDs were calculated between the 1D-FDMs and the 1D Fourier solutions. Fig. S6c shows the 1D-FDMs MSD 21 has initial MSD peaks comparable to the 1D-uFDMs initial MSD peaks for ∆𝑦 = 0.02, 0.04 𝜇𝑚 , with the initial MSD peaks increasing in size with ∆𝑥 (compare Fig. S6a, S6b and S6c). As ∆𝑦 increases further interpolation inaccuracies are combined with the 1D-FDMs inherent error to increase the initial 1D-uFDMs MSD peaks, nonetheless these increased peak sizes remain within the same order of magnitude as the 1D-FDMs MSD peaks. To further confirm the hypothesis that the 1D-uFDMs initial MSD peaks could be attributed to the 1D-FDMs inherent error we compared 1D-FDMs and 1D-uFDMs SD kymographs of the first 20 s of SD data (Fig. S6d). Indeed, the 1D-uFDMs SD kymographs showed close resemblance to the 1D-uFDM SD kymographs for smaller ∆𝑦 , with the largest SDs at the points of greatest curvature in the concentration profile. The last point of note is that the 1D-FDMs MSDs do not increase again after the initial peak (Fig. S6c) confirming that the dominant error increase in 1D-uFDMs MSD is due to interpolation inaccuracy at the domain edges (Fig. S5). Simulation set 2: In the steady state analysis it was found that increasing ∆𝑡 improved the accuracy of the semi-implicit 1D-uFDM when estimating the homogeneous steady state, but not the explicit 1D-uFDM (Fig. S4b). To ask if the error inherent in the 1D-FDMs can account for the differences between the explicit and semi-implicit 1D-uFDM MSDs at steady state we performed a time course MSD analysis on 1D-FDMs and 1D-uFDMs with fixed ∆𝑥 = ∆𝑦 = 0.1 𝜇𝑚 and increased ∆𝑡 from (the ∆𝑡 value in simulation set 1) to for explicit 1D-FDM and 1D-uFDM and to for the implicit 1D-FDM and semi-implicit 1D-uFDM, ensuring a range of values for ∆𝑡 while retaining numerical stability. We found that the error inherent in 1D-FDMs followed the same trend as did the steady state 1D-uFDMs MSD, namely increasing ∆𝑡 in the explicit 1D-FDM had no effect on the MSD, and increasing ∆𝑡 in the implicit 1D-FDM had a dramatic effect on the MSD (compare Fig. S4b and S6e). The 1D-FDM and 1D-uFDM initial MSD peaks followed the same trend, as expected (compare Fig. S6e and S6f inset). Looking at the entire MSD time course for the semi-implicit 1D-uFDM (Fig. S6f) we see that, for the values of ∆𝑡 chosen, the increase in MSD initial peak offsets the steady state MSD. Simulation set 3: For completion we analysed the MSD dynamics in the semi-implicit 1D-uFDM for fixed Δ𝑥 = 0.1 𝜇𝑚 , increasing
Δ𝑦 = 0.01, 0.02, … 0.5 𝜇𝑚 , and increasing ∆𝑡 − (Δ𝑦, 𝐷, 𝛽) =∆𝑡 − (Δ𝑦, 0.1, 0.1) = ∆𝑦 𝑠 (steady state MSD shown in Fig. S4f). Increasing Δ𝑦 (and Δ𝑡 ) resulted in an increasing MSD initial peak, it did not have a minimum similar to the MSD at steady state, indicating the optimal Δ𝑦 for minimum MSD at steady state is found at a balance point between 1D-FDM error and interpolation boundary error (Fig. S6g). Appendix 8: Accuracy of the semi-implicit 1D-uFDM when simulating diffusion.
Fig. 2B-E shows the comparisons between the explicit 2D-FDM, 1D-FDM and 1D-uFDM when modelling diffusion in the full and reduced-dimension models. Fig. S7 shows the same investigation, but this time comparting the results of the explicit 2D-FDM, 1D-FDM and semi-implicit 1D-uFDM.
𝐷 = 0.1 𝜇𝑚 𝑠 −1 . A set of Δ𝑥 , Δ𝑦 and Δ𝑡 values were chosen that gave an ‘improved’ steady state accuracy with a fairly small initial MSD peak, namely Δ𝑥 = Δ𝑦 = 0.1 𝜇𝑚 , Δ𝑡 = 0.1 𝑠 (Fig. S6f). As expected setting
Δ𝑡 = 0.01 𝑠 , the same as Δ𝑡 in the explicit comparison, gave a semi-implicit curve comparable to the explicit MSD curve (compare Fig. S7b with Fig. 2E). Appendix 9: Accuracy of the 1D-uFDM and 1D-uFDMs when simulating FRAP.
The accuracy analysis for the 1D diffusion solutions when solving for an initial peak or an initial trough (i.e. the FRAP solution) is almost identical (compare Fig. 2C and E with Fig. S8b and c). Both the inherent 1D-FDM error and edge error (Appendix 7) increase proportionally with the concentration’s curvature, and in both the peak and trough profiles the curvature is described by the same equation, 𝑒 −(𝑥 +𝑦 ) in the 2D case, 𝑒 −𝑥 for 1D. Both peak and trough solutions are 22 subject to the same diffusion coefficient so the curvature of the profiles change at the same rate giving the same accuracy analysis results. Appendix 10: Reaction-diffusion comparisons and the effect of geometry.
The explicit 1D-uFDM captures RD dynamics through the focal plane more accurately than the explicit 1D-FDM (Fig. S9), with an MSD an order of magnitude smaller than that of the 1D-FDM for all 𝛼 (Fig. S9c and f). Analysis of implicit 1D-FDM and semi-implicit 1D-uFDM comparisons gave almost identical results to the explicit comparisons, data not shown. Recall, if the zero flux assumption holds, then the 1D-FDM solution is an accurate approximation of molecular dynamics through the focal plane of the 2D membrane. Thus the 1D-FDM solution to the RD equations is accurate if we were investigating an initial pulse on the surface of a spherical cell, for example (Fig. 1A, E). Fig. S10 shows the RD equations, 𝛼 = 3 , solved on the body of an elongated cell using the square mesh and on a spherical cell using a spherical mesh. The same RD equations produce a patch on the square mesh and a ring on the spherical mesh. 23 Figures and Tables
Fig. S1: Concentration profiles used to investigate the accuracy of estimating concentrations at phantom points. (a)
1D concentration profiles described by the equation, 𝒖 =𝒆 −𝜶𝒙 . (b) Example of the 2D system before dimension reduction reduced, 𝜶 = 𝟏 , 𝒖 = 𝒆 −(𝒙 +𝒚 ) . (c)
1D concentration profiles described by the equation, 𝒖 = 𝟏 − 𝒆 −𝜶𝒙 . (d) Example of the 2D system before dimension reduction, 𝜶 = 𝟏 , 𝒖 = 𝟏 − 𝒆 −(𝒙 +𝒚 ) . (e)
1D concentration profiles described by the equation, 𝒖 = 𝐦𝐢𝐧(𝒆 −𝜶 𝒙 , 𝟏 − 𝒆 −𝜶 𝒙 ) . 𝜶 = 𝟏𝟎 for all curves. 𝜶 = 𝟏, 𝟎. 𝟏, 𝟎. 𝟎𝟏 for the yellow, red, blue curve, respectfully. (f) Example of the 2D system before dimension reduction, 𝜶 = 𝟏 , 𝒖 = 𝐦𝐢𝐧(𝒆 −(𝒙 +𝒚 ) , 𝟏 − 𝒆 −𝟏𝟎(𝒙 +𝒚 ) ) . a c b e d f -6 space ( 𝜇𝑚 ) s p a c e ( 𝜇 𝑚 ) c o n c e n t r a t i o n ( 𝜇 𝑀 ) -6 space ( 𝜇𝑚 ) s p a c e ( 𝜇 𝑚 ) c o n c e n t r a t i o n ( 𝜇 𝑀 ) -6 space ( 𝜇𝑚 ) s p a c e ( 𝜇 𝑚 ) c o n c e n t r a t i o n ( 𝜇 𝑀 ) Fig. S2: SSD between estimated concentrations at phantom points and actual concentrations in the 2D System.
The analysis was performed on the profiles, 𝒖 = 𝒆 −𝜶𝒙 and 𝒖 = 𝒆 −𝜶(𝒙 +𝒚 ) (Fig. S1a and b). 25 Fig. S3: Numerical stability boundaries.
Black lines show the numerical stability boundaries for ∆𝒚 = 𝟎. 𝟎𝟏, 𝟎. 𝟎𝟐, … , 𝟎. 𝟏 𝝁𝒎 . 𝑫 = 𝟎. 𝟏 𝝁𝒎 𝒔 −𝟏 Values of ∆𝒕 on or below the black lines are numerically stable. Green dots show the stability boundary when ∆𝒙 = ∆𝒚 . Note, the ∆𝒕 axis have different maxima on each graph. (a) Boundaries for numerical stability in the explicit 1D-uFDM are given by the equation, ∆𝒕 = ( ∆𝒙 ∆𝒚 +∆𝒙 ) 𝒔 . (b) Boundaries for numerical stability in the semi-implicit 1D-uFDM are given by the equation, ∆𝒕 = ∆𝒚 𝑫 𝒔 . 26 Fig. S4: Accuracy of unrestricted movement finite difference methods at steady state.
Solid lines show the results of the explicit method, dotted the semi-implicit. In (a), (c) and (e) Lighter grey indicates increasing ∆𝒙 or ∆𝒚 . (a) MSD between the explicit and semi-implicit 1D-FDMs and the 2D 27 Fourier solutions at steady state, for ∆𝒙 = 𝟎. 𝟎𝟏, 𝟎. 𝟎𝟐, … , 𝟎. 𝟏 𝝁𝒎 , ∆𝒚 = 𝟎. 𝟎𝟏, 𝟎. 𝟎𝟐, … , 𝟎. 𝟏 𝝁𝒎 and 𝚫𝒕 = 𝟎. 𝟎𝟎𝟎𝟑 . (b) MSD at steady state for ∆𝒙 = 𝚫𝒚 = 𝟎. 𝟏 𝝁𝒎 , varying ∆𝒕 . (c) MSD between the explicit and semi-implicit 1D-FDMs and the 2D Fourier solutions at steady state, for ∆𝒙 ≥𝟎. 𝟎𝟓, 𝚫𝒚 ≥ 𝟎. 𝟎𝟓 and ∆𝒕 = 𝟎. 𝟎𝟎𝟕𝟓 . (d) MSD at steady state for ∆𝒙 = 𝚫𝒚 = 𝟎. 𝟓 𝝁𝒎 , varying ∆𝒕 . (e) MSD at steady state for ∆𝒙 = 𝟎. 𝟎𝟏, 𝟎. 𝟎𝟐, … , 𝟎. 𝟏 𝝁𝒎, ∆𝒚 = 𝟎. 𝟎𝟏, 𝟎. 𝟎𝟐, … , 𝟎. 𝟏 𝝁𝒎 , and ∆𝒕 − = ( ∆𝒙 ∆𝒚 +∆𝒙 ) 𝒔 for explicit simulations and ∆𝒕 − = ∆𝒚 𝒔 for semi-implicit. (f) MSD at steady state for semi-implicit 1D-uFDM, ∆𝒙 = 𝟎. 𝟏 𝝁𝒎 , ∆𝒚 = 𝟎. 𝟎𝟏, 𝟎. 𝟎𝟐, … 𝟎. 𝟓 𝝁𝒎 , and ∆𝒕 − (𝚫𝒚, 𝟎. 𝟏, 𝟎. 𝟏) = ∆𝒚 𝒔 . 28 Fig. S5: Accuracy dynamics of unrestricted movement finite difference methods: Interpolation error. (a)
MSD between the explicit and semi-implicit 1D-FDMs and the 2D Fourier solutions for ∆𝒙 = 𝟎. 𝟎𝟐, 𝟎. 𝟎𝟒, … , 𝟎. 𝟏 𝝁𝒎 , ∆𝒚 = 𝟎. 𝟎𝟐, 𝟎. 𝟎𝟒, … , 𝟎. 𝟏 𝝁𝒎 and 𝚫𝒕 = 𝟎. 𝟎𝟎𝟎𝟑 𝒔 . For each 𝚫𝒚 the five MSD curves for ∆𝒙 = 𝟎. 𝟎𝟐, 𝟎. 𝟎𝟒, … , 𝟎. 𝟏 𝝁𝒎 are drawn on top of each other. (b) Cartoon showing the estimation of 𝒖 √(𝑵/𝟐 ∆𝒙) +∆𝒚 𝝉 at the boundary. (c) Kymograph of the absolute distance between each point of the semi-implicit 1D-uFDM and the central row of the 2D Fourier solution, at every second, for
𝚫𝒙 = 𝟎. 𝟏 𝝁𝒎 , 𝚫𝒚 = 𝟎. 𝟎𝟐 𝝁𝒎 , 𝚫𝒕 = 𝟎. 𝟎𝟎𝟎𝟑 𝒔 . (d) Kymograph, with upper threshold defined by the steady state concentration, of the concentration profile of the semi-implicit 1D-uFDM for
𝚫𝒙 = 𝟎. 𝟏 𝝁𝒎 , 𝚫𝒚 = 𝟎. 𝟎𝟐 𝝁𝒎 , 𝚫𝒕 = 𝟎. 𝟎𝟎𝟎𝟑 𝒔 . (e) Graph of the absolute distance and concentration profile at the boundary. 29 30
Fig. S6: Accuracy dynamics of unrestricted movement finite difference methods: 1D-FDM error. (a)
Initial MSD peaks between the explicit 1D-uFDM and the 2D Fourier solutions, for ∆𝒙 =𝟎. 𝟎𝟐, 𝟎. 𝟎𝟒, … , 𝟎. 𝟏 𝝁𝒎 , ∆𝒚 = 𝟎. 𝟎𝟐, 𝟎. 𝟎𝟒, … , 𝟎. 𝟏 𝝁𝒎 and 𝚫𝒕 = 𝟎. 𝟎𝟎𝟎𝟑 𝒔 . (b) Values of initial MSD peaks for explicit and semi-implicit 1D-uFDM, for ∆𝒙 = 𝟎. 𝟎𝟐, 𝟎. 𝟎𝟒, … , 𝟎. 𝟏 𝝁𝒎 , ∆𝒚 =𝟎. 𝟎𝟐, 𝟎. 𝟎𝟒, … , 𝟎. 𝟏 𝝁𝒎 and 𝚫𝒕 = 𝟎. 𝟎𝟎𝟎𝟑 𝒔 . (c) MSD between the explicit 1D-FDM and the 1D Fourier solutions, for ∆𝒙 = 𝟎. 𝟎𝟐, 𝟎. 𝟎𝟒, … , 𝟎. 𝟏 𝝁𝒎 and
𝚫𝒕 = 𝟎. 𝟎𝟎𝟎𝟑 𝒔 . Lighter grey indicates larger 𝚫𝒙 , key in Fig. S6a. (c inset) Values of initial MSD peaks for explicit and semi-implicit 1D-FDM. (d)
SD kymographs showing the initial SD peak for ∆𝒙 = 𝟎. 𝟎𝟐, 𝟎. 𝟎𝟒, … , 𝟎. 𝟏 𝝁𝒎 , 𝚫𝒕 = 𝟎. 𝟎𝟎𝟎𝟑 𝒔 . Row 1 explicit 1D-FDM, rows 2 to 6 explicit 1D-uFDM for ∆𝒚 = 𝟎. 𝟎𝟐, 𝟎. 𝟎𝟒, … , 𝟎. 𝟏 𝝁𝒎 . (e) Initial MSD peaks for 1D-FDM, ∆𝒙 = 𝚫𝒚 = 𝟎. 𝟏 𝝁𝒎 , varying ∆𝒕 . (f) MSD for semi-implicit 1D-uFDM, ∆𝒙 = 𝚫𝒚 =𝟎. 𝟏 𝝁𝒎 , varying ∆𝒕 . (f inset) Initial MSD peaks for 1D-uFDM, ∆𝒙 = 𝚫𝒚 = 𝟎. 𝟏 𝝁𝒎 , varying ∆𝒕 . (g) MSD for semi-implicit 1D-uFDM, ∆𝒙 = 𝟎. 𝟏 𝝁𝒎 , ∆𝒚 = 𝟎. 𝟎𝟏, 𝟎. 𝟎𝟐, … 𝟎. 𝟓 𝝁𝒎 , and ∆𝒕 − (𝚫𝒚, 𝟎. 𝟏, 𝟎. 𝟏) = ∆𝒚 𝒔 . (g inset) Initial MSD peaks for semi-implicit 1D-uFDM, ∆𝒙 = 𝟎. 𝟏 𝝁𝒎 , ∆𝒚 = 𝟎. 𝟎𝟏, 𝟎. 𝟎𝟐, … 𝟎. 𝟓 𝝁𝒎 , and ∆𝒕 − (𝚫𝒚, 𝟎. 𝟏, 𝟎. 𝟏) = ∆𝒚 𝒔 . 31 Fig. S7: Accuracy dynamics of the implicit 1D-FDM and semi-implicit 1D-uFDMs when simulating diffusion in reduced dimension models. (a)
Kymographs of molecule 𝒖 on the focal plane of the 2D model and in the 1D models. (b) Maximum concentration dynamics. (c)
Mean concentration in the 1D solutions compared with the mean concentration in the focal plane of the 2D solution. (d)
MSD between the focal plane in the 2D solution and the 1D solutions. 1D-FDM comparison shown in blue, 1D-uFDM comparison shown in red. 32
Fig. S8: Accuracy dynamics of the 1D-FDM and 1D-uFDMs when simulating FRAP in reduced dimension models.
Figures support Fig. 3B-E. (a)
Time course of the minimum concentration, i.e. the concentration at the centre of the ROI. (b)
The MSD between the explicit 2D solution and the explicit 1D-FDM, explicit 1D-uFDM and semi-implicit 1D-uFDM. (c)
Kymographs of the SD between the explicit 2D solution and the explicit 1D-FDM, explicit 1D-uFDM and semi-implicit 1D-uFDM. 33
Fig. S9: Accuracy dynamics of the 1D-FDM and 1D-uFDMs when solving RD equations.
Figures support Fig. 4B-D. (a)
Mean concentration in the 1D solutions compared with the mean concentration in the focal plane of the 2D solution. (b) 𝒖 steady state colorplots for the 2D solution on the focal plane, and 1D solutions. 𝜶 = 𝟎. 𝟓 in panels (c), (d) and (e). (c) Kymographs of the 𝒖
34 and 𝒖 concentrations in the focal plane of the explicit 2D solution and 1D solutions. (d) 𝒖 and 𝒖 concentration profiles at steady state. (e) 𝒖 and 𝒖 MSD dynamics, comparing the concentrations in the focal plane of the 2D solution and the 1D solutions. 𝜶 = 𝟑 in panels (f), (g) and (h). (f)
Kymographs of the 𝒖 and 𝒖 concentrations in the focal plane of the 2D solution and 1D solutions. (g) 𝒖 and 𝒖 concentration profiles at steady state. (h) 𝒖 and 𝒖 MSD dynamics, comparing the concentrations in the focal plane of the 2D solution and the 1D solutions. Fig. S10: The effect of geometry on the RD solutions. Compare with Fig. 4C-D , 𝜶 = 𝟑 . (a)
Final steady state of the 2D-FDM solution to the RD equations, initial condition 𝜶 = 𝟑 , solved on the body of an elongated cell. (b)
Final steady state of the 2D-FDM solution to the RD equations, initial condition 𝜶 = 𝟑 , solved on a spherical mesh using finite differences. a. b. 𝑢 c o n c e n t r a t i o n ( µ M ) Table S1: Minimum and maximum SSDs between estimated concentrations at phantom points and actual concentrations in the 2D System.
Profiles analyzed are those in Fig. S1a, c and e. Minimum SSDs shown in blue text, maximum in black. Table S2: Results of the numerical test of the explicit 1D-uFDM stability condition.
Smallest values of 𝛃 out of those tested ( 𝛃 = 𝟏𝐱𝟏𝟎 −𝟒 , 𝟏𝐱𝟏𝟎 −𝟑 , 𝟏𝐱𝟏𝟎 −𝟐 ) for which simulations with ∆𝒕 + were unstable and ∆𝒕 − were stable. 𝑫 = 𝟎. 𝟏 𝝁𝒎 𝒔 −𝟏 . Blue cells indicate the values of ∆𝒙, ∆𝒚 for which ∆𝒕 − simulations with 𝜷 = 𝟏𝐱𝟏𝟎 −𝟒 , 𝟏𝐱𝟏𝟎 −𝟑 were numerically unstable. ∆𝒙 ∆𝒚 Table S3: Results of the numerical test of the explicit 1D-uFDM stability condition. ∆𝒕 ± =𝚫𝒕 ± 𝜷𝚫𝒕 for 𝜷 = 𝟎. 𝟎𝟏 . ∆𝒙 ∆𝒚 -4 ± 4.7x10 -6 -4 ± 4.8x10 -6 -4 ± 4.9x10 -6 -4 ± 4.9x10 -6 -3 ± 1.9x10 -5 -4 ± 4.9x10 -6 -3 ± 1.9x10 -5 -4 ± 5.0x10 -6 -3 ± 1.9x10 -5 -4 ± 5.0x10 -6 -3 ± 2.0x10 -5 -3 ± 4.2x10 -5 -4 ± 5.0x10 -6 -3 ± 2.0x10 -5 -3 ± 4.3x10 -5 Table S4: Results of the numerical test of the semi-implicit 1D-uFDM stability condition.
Smallest values of 𝛃 out of those tested ( 𝛃 = 𝟏𝐱𝟏𝟎 −𝟒 , 𝟏𝐱𝟏𝟎 −𝟑 , 𝟏𝐱𝟏𝟎 −𝟐 ) for which simulations with ∆𝒕 + were unstable and ∆𝒕 − were stable. 𝑫 = 𝟎. 𝟏 𝝁𝒎 𝒔 −𝟏 . ∆𝒙 ∆𝒚 Table S5: Summary of the 1D-uFDM accuracy analysis.
Full analysis in Appendix 6 and 7. ↑ 𝚫𝒚 slight increase ↑ 𝚫𝒙 slight reduction ↑ 𝚫𝒕 no difference reduction ↑ 𝚫𝒚 & 𝚫𝒕 slight increase great reduction, up to a point
Initial MSD Peak ↑ 𝚫𝒚 slight increase ↑ 𝚫𝒙 slight increase for
Δ𝑦 = 0.02 to 0.06 slight reduction for
Δ𝑦 = 0.1 ↑ 𝚫𝒕 no difference increase ↑ 𝚫𝒚 & 𝚫𝒕
N/A great increase 41
References
1. A. M. Turing, The Chemical Basis of Morphogenesis.
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences , 37-72 (1952). 2. L. Wolpert, Positional Information and the Spatial Pattern of Cellular Differentiation.
Journal of Theoretical Biology , 1-47 (1969). 3. R. J. Payne, C. S. Grierson, A theoretical model for ROP localisation by auxin in Arabidopsis root hair cells. PLoS One , e8337 (2009). 4. N. W. Goehring et al. , Polarization of PAR proteins by advective triggering of a pattern-forming system. Science , 1137-1141 (2011). 5. S. E. Smith et al. , Independence of symmetry breaking on Bem1-mediated autocatalytic activation of Cdc42.
J Cell Biol , 1091-1106 (2013). 6. V. Breña-Medina, A. R. Champneys, C. Grierson, M. J. Ward, Mathematical Modeling of Plant Root Hair Initiation: Dynamics of Localized Patches.
SIAM Journal on Applied Dynamical Systems , 210-248 (2014). 7. S. Blanchoud, C. Busso, F. Naef, P. Gonczy, Quantitative analysis and modeling probe polarity establishment in C. elegans embryos. Biophys J , 799-809 (2015). 8. J. G. Chiou et al. , Principles that govern competition or co-existence in Rho-GTPase driven polarization.
PLoS Comput Biol , e1006095 (2018). 9. M. Onsum, C. V. Rao, A mathematical model for neutrophil gradient sensing and polarization. PLoS Comput Biol , e36 (2007). 10. R. H. William, L. Benjamin, L. Andre, E.-K. Leah, Modelling cell polarization driven by synthetic spatially graded Rac activation. PLoS Computational Biology , e1002366-e1002366 (2012). 11. P. Rangamani et al. , Decoding Information in Cell Shape. Cell , 1356-1369 (2013). 12. N. Muller et al. , A Predictive Model for Yeast Cell Polarization in Pheromone Gradients.
PLoS Comput Biol , e1004795 (2016). 13. R. Gessele, J. Halatek, L. Wurthner, E. Frey, Geometric cues stabilise long-axis polarisation of PAR protein patterns in C. elegans. Nat Commun , 539 (2020). 14. H. Meinhardt, Tailoring and coupling of reaction-diffusion systems to obtain reproducible complex pattern formation during development of the higher organisms. Applied Mathematics and Computation , 103 (1989). 15. A. T. Dawes, L. Edelstein-Keshet, Phosphoinositides and Rho Proteins Spatially Regulate Actin Polymerization to Initiate and Maintain Directed Movement in a One-Dimensional Model of a Motile Cell. Biophysical Journal , 744-768 (2007). 16. Y. Wu et al. , Rapid diffusion-state switching underlies stable cytoplasmic gradients in the Caenorhabditis elegans zygote. Proc Natl Acad Sci U S A , E8440-E8449 (2018). 17. L. V. Fausett,
Applied numerical analysis using MATLAB (Pearson Prentice Hall, ed. 2nd ed., 2008). 18. M. Kang, C. A. Day, A. K. Kenworthy, E. DiBenedetto, Simplified equation to extract diffusion coefficients from confocal FRAP data.
Traffic , 1589-1600 (2012). 19. J. Crank, The mathematics of diffusion / by J. Crank , Oxford science publications. (Clarendon Press, Oxford [England], 1975). 20. S. A. Ramirez, D. J. Lew, S. Raghavachari, Dendritic spine geometry can localize GTPase signaling in neurons.
Molecular Biology of the Cell , 4171-4181 (2015).
21. C. T. Lee et al. , 3D mesh processing using GAMer 2 to enable reaction-diffusion simulations in realistic cellular geometries.
PLoS Comput Biol , e1007756 (2020). 22. W. E. Boyce, R. C. DiPrima, Elementary differential equations and boundary value problems (J. Wiley, ed. 6th ed., 1996). 23. K. A. Stroud,
Fourier series and harmonic analysis (Thornes, 1984).(Thornes, 1984).