Describing the movement of molecules in reduced-dimension models

When addressing spatial biological questions using mathematical models, symmetries within the system are often exploited to simplify the problem by reducing its physical dimension. In a reduced-dimension model molecular movement is restricted to the reduced dimension, changing the nature of molecular movement. This change in molecular movement can lead to quantitatively and even qualitatively different results in the full and reduced systems. Within this manuscript we discuss the condition under which restricted molecular movement in reduced-dimension models accurately approximates molecular movement in the full system. For those systems which do not satisfy the condition, we present a general method for approximating unrestricted molecular movement in reduced-dimension models. We will derive a mathematically robust, finite difference method for solving the 2D diffusion equation within a 1D reduced-dimension model. The methods described here can be used to improve the accuracy of many reduced-dimension models while retaining benefits of system simplification. Natasha Savage develops a general finite difference method for modeling 2D diffusion using a 1D reduced model. This method can be useful to improve the accuracy of many reduced-dimension models while retaining benefits of system simplification.

S imple 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][4][5][6][7][8] (Fig. 1a). When addressing questions which include cytoplasmic gradients, for example, a system's dimension is often reduced from 3D to 2D 1,[9][10][11][12][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 and c). Thus, when one reduces the dimension of a system and restricts molecular movement to the reduced dimension, they are changing the geometry of the problem. It is understood that cell geometry influences molecular movement and patterning 11,13,[17][18][19] . Here we present a general methodology that enables a reduced-dimension model to take the system's full geometry into account, by using it to estimate the movement of molecules through the focal plane. We go on to use this general methodology to derive a simple numerical method (the 1D-uFDM) that can be used to solve the 2D diffusion equation in a 1D reduced-dimension model.

Results
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. Throughout the manuscript 2D and 1D solutions are compared to analyse and illustrate the accuracy of the reduced-dimension method being proposed. A user of this method would not generate the full dimension model, they would only generate the simplified reduced-dimension model.
Reduced-dimension models have a zero-flux assumption. In order to use a numerical method to calculate changes in molecular concentrations over space and time, the space and time must be discretised. For a 2D space, discretisation is achieved by drawing a mesh over the space (Fig. 1a, d, e). The concentration within each discretised space is represented by a single value and this value is used in the numerical method. Discretisation of time is analogous to collecting time course data, one defines a time step, 4t seconds, and collects data at a number of time points 4t seconds apart. The finite difference method is a commonly used method for calculating molecular concentrations over discretised space and time 20,21 .
To investigate any inbuilt assumptions of reduced-dimension models we compared the finite difference solution to the diffusion equation in 2D and 1D. For the 2D solution cover the 2D membrane with a regular mesh (Fig. 1d). The mesh points along the x-axis are labelled i ¼ 1; 2; ; N and are distance 4x μm apart, the mesh points along the y-axis are labelled j ¼ 1; 2; ; M and are distance 4y μm apart. Let u τ i;j represent the concentration of molecule u on the membrane at point i; j À Á and time point τ, τ ¼ 0; 1; 2; 3; ; time points are 4t seconds apart. The change in concentration of u τ i;j 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 D μm 2 s À1 , is 20,21 Explicit 2D-FDM: The finite difference method is an iterative method. Thus, one must know the concentrations of molecule u at all points of the Fig. 1 Dimension reduction. a Dimension reduction example showing a polarity patch on a spherical cell and body of an elongated cell with their reduceddimension 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. The focal plane is shown as a white stripe. c Molecules in a reduced-dimension model are restricted to the focal plane (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, u τ i;JÀ1 ¼ u τ i;J ¼ u τ i;Jþ1 . f Regular mesh over a polarity patch on the body of an elongated cell. Distances from the centre of the patch of each point in rows J ± 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 e u τ i;J ± 1 (black open circles). mesh at time point τ ¼ 0. We denote these initial concentrations as u 0 i;j . The initial concentrations (u 0 i;j ) are substituted into Eq. (1) to calculate u 1 i;j . u 1 i;j are the concentrations at time point 1. The time at time point 1 is Δt seconds. The concentrations u 1 i;j are then substituted into Eq. (1) to calculate u 2 i;j , the concentrations at time point 2, 2Δt seconds. This iterative process continues until adequate time course data are calculated. The calculated time course data are referred to as the solution to the finite difference method.
Note that in order to calculate the concentrations u τþ1 i;j the 2D-FDM calculates molecular movement along the x-axis and y-axis separately. Molecular movement along the x-axis is calculated by the term 4t , molecular movement along the y-axis is calculated by the term 4t Á . To build a 1D reduced-dimension model describing the diffusive movement of u on a 2D surface, the focal plane is set to run through the axis of symmetry of u 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 x-axis at row j ¼ J (Fig. 1d). The explicit 1D finite difference scheme (1D-FDM) used to solve the 1D diffusion equation is 20,21 Explicit 1D-FDM: As the 1D reduced-dimension model only contains information about concentrations on row J the 1D-FDM contains no terms for calculating the movement of u along the y-axis (compare the 2D and 1D-FDMs, Eqs. 1 and 2). 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 Á ¼ 0, the number of molecules leaving the focal plane, À2u τ i;J , is equal to the number of molecules entering it, u τ i;JÀ1 þ u τ i;Jþ1 , at all points in space, i, for all time, τ: There is zero-flux through the focal plane.
Conditions under which the zero-flux assumption is valid. The zero-flux assumption is valid for a subset of reduced-dimension models, those for which a mesh can be drawn such that the zeroflux 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 u τ i;JÀ1 , u τ i;J and u τ i;Jþ1 are equal and the zero-flux assumption holds (Fig. 1e). Thus, molecular movement in a 1D reduced-dimension model of this system would be accurately described by the 1D diffusion equation.
An example of a system for which the zero-flux assumption does not hold is 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). Thus, the 1D diffusion equation would give an inaccurate approximation of 2D molecular movement through the focal plane.
For the modeller, the validity of the zero-flux assumption can be ascertained without generating a solution for the full system and calculating molecular movement through the focal plane. The modeller can consider the symmetries in the full system and the form of the full mesh using cartoons, as in Fig. 1. From the cartoons one can estimate whether or not the zero-flux assumption holds. This cartoon estimation of zero-flux is not dissimilar to the estimation made by the modeller when deciding whether or not it is appropriate to reduce a system's dimension.
A general method for calculating molecular movement into and out of the focal plane. Here we present a general methodology which can be used to increase the accuracy of reduceddimension models that do not satisfy the zero-flux assumption: If it is possible to estimate the concentrations either side of the focal plane (for example, terms u τ i;Jþ1 and u τ i;JÀ1 in the explicit 2D-FDM equation) using the concentrations on the focal plane (u τ i;J in our example), then we can estimate molecular movement into and out of the focal plane. In order to reduce a system's dimension, the system must exhibit symmetry. Therefore, by definition, the concentration profile in the focal plane (i.e. in the dimensionally reduced model) represents the concentration profile of the full system and thus holds information about the concentrations at every point in the full system. As such, the information in the focal plane can be used to calculate concentrations either side of the focal plane in any system for which dimension reduction is possible.
Below we provide a proof of principle for the general method. We derive a mathematically robust numerical method which can be used to improve the accuracy of a subset of reduced-dimension models. Namely, diffusive movement in 1D reduced-dimension models that do not satisfy the zero-flux assumption and exhibit radial molecular dynamics. This proof of principle does not represent the limits of the general method for estimating molecular movement into and out of the focal plane. As discussed above, the general method can be used for systems of any dimension and concentration profiles that do not exhibit radial symmetry. Furthermore, molecular movement need not be governed by diffusion.
Solving the 2D diffusion equation in a 1D reduced-dimension model. We will use the general methodology to construct a FDM that solves the 2D diffusion equation in a 1D reduced-dimension model. Consider again the polarity patch on the body of an elongated cell (Fig. 1a, d, f). The 2D polarity patch has radial symmetry. A property of radial symmetry is that the concentration profile along all lines running through the centre of the patch is identical. When we reduce the dimension of this system to 1D, we position the focal plane such that the 1D model calculates the concentration profile along one line running through the centre of the 2D patch (row J). As the concentration profile on all lines running through the centre of the patch is identical, we can use the concentration profile in the 1D model to calculate the concentrations at all points on the 2D surface. The explicit 2D-FDM (Eq. 1) contains two concentrations that are not on the focal plane, u τ i;Jþ1 and u τ i;JÀ1 . Thus, to calculate molecular movement into and out of the focal plane, row J, and solve the 2D diffusion equation in a 1D reduced-dimension model, we need only estimate concentrations on either side of the focal plane, rows J ± 1.
Letũ τ i;J ± 1 denote the estimated concentrations on rows J ± 1. u τ i;J ± 1 can be estimated from the concentrations of molecule u along the focal plane, u τ i;J , using interpolation. The interpolation mesh points are found using Pythagoras' theorem ( Fig. 1f, g, "Methods: Generating the interpolation mesh in 1D"). In order to use Pythagoras to calculate the interpolation mesh points, the modeller has to set a value for Δy. In the illustrative examples below, we have set Δy ¼ Δx. Analysis on the accuracy of estimating concentrations either side of the focal plane and the choice of Δx, Δy, can be found in Supplementary: Accuracy of estimating concentrations at phantom points using interpolation.
Further discussion on the derivation of the explicit 1D finite difference diffusion equation with unrestricted movement (explicit 1D-uFDM) can be found in "Methods: Explicit 1D-uFDM derivation". The 1D-uFDM is Explicit 1D-uFDM: See "Methods: Explicit 1D-uFDM solution" for the solution to the explicit 1D-uFDM. The explicit 1D-uFDM numerically stability condition is derived in "Methods: Explicit 1D-uFDM numerical stability condition" and tested numerically in Supplementary: Numerical verification of stability Conditions. A fully implicit 1D-uFDM is ill defined (Supplementary: An implicit 1D-uFDM is ill defined); however, a semi-implicit 1D-uFDM and numerical stability condition can be derived ("Methods: Semiimplicit 1D-uFDM solution, Methods: Semi-implicit 1D-uFDM numerical stability condition"). The explicit and semi-implicit 1D-uFDMs solve the 2D diffusion equation in a 1D reduceddimension model.
In all three illustrative examples below, we are testing the ability of the explicit 1D-uFDM to solve the 2D diffusion equation in a 1D reduced-dimension model. We will also illustrate the accuracy gained by solving the 2D diffusion equation in a 1D reduced-dimension model using the 1D-uFDM, when compared to using the 1D diffusion equation. Thus, the solutions of the 1D-uFDM and 1D-FDM will be compared to the solution on a slice through the centre of the patch in the full system. The 2D solution will be calculated using 2D-FDM.
Illustrative example 1: Diffusion. First, we considered the diffusion of molecules, u, from a concentrated patch on the membrane. The initial concentration profile in the reduced-dimension models (Fig. 2a) was identical to the initial concentration along a slice through the centre of the patch in the 2D system (Fig. 2a,  inset). The diffusion equation was solved using the explicit 2D-FDM, 1D-FDM and 1D-uFDM (Eqs. 1-3), diffusion coefficient D ¼ 0:1 μm 2 s À1 ("Methods: Illustrative example 1: Parameters). The solutions of u in both 1D reduced-dimension models were compared with the solutions of u along a slice through the centre of the initial patch (see Supplementary: Accuracy of the semiimplicit 1D-uFDM when simulating diffusion 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 (1D-FDM, Eq. 2) was quantitatively inaccurate at estimating the concentration on a slice through the centre of the 2D patch (Fig. 2b, c). This is because on the 2D membrane molecules diffuse out of the focal plane resulting in a reduction in the mean concentration of molecules in the central slice (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 1D-uFDM estimates the movement of molecules out of the focal plane producing a more accurate reduced-dimension representation of the full system (Fig. 2b-e). While the 1D-uFDM represents the 2D system well it does contain error (Fig. 2e), in-depth error analysis can be found in Supplementary: Steady state accuracy of 1D-uFDMs and Accuracy dynamics of 1D-uFDMs.
Illustrative example 2: Florescence recovery after photobleaching. The 1D-uFDM was able to solve the 2D diffusion equation in a 1D reduced-dimension model when calculating molecules diffusing away from a central patch. We went on to ask, how accurate is the 1D-uFDM when calculating the movement of molecules into a central trough. To answer this question we modelled a florescence recovery after photobleaching (FRAP) experiment. FRAP experiments are used to estimate the diffusion coefficient of a fluorescently tagged protein. In a FRAP experiment fluorophores, attached to a protein of interest, are bleached by a laser. The florescence recovery within the bleached area is recorded and used to estimate a diffusion coefficient. In our FRAP model the initial fluorophore concentrations were set to reflect a uniformly covered membrane after bleaching with a Gaussian laser 22 (Fig. 3a). Diffusive movement of the unbleached, fluorescently tagged, proteins was modelled using the 2D-FDM, the 1D-FDM and the 1D-uFDM (Eqs. 1-3), diffusion coefficient D ¼ 0:1 μm 2 s À1 . Results: In all FRAP solutions fluorescently tagged proteins diffused into the bleached area (Fig. 3b). The 1D-FDM solution had a low homogeneous steady-state concentration when compared to the homogeneous steady state along a slice through the centre of the bleached area in the full system (Fig. 3b). The low steady-state concentration in the 1D-FDM solution can be accounted for by the inability of tagged proteins to move into the focal plane (Fig. 3c). The 1D-uFDM estimated the movement of tagged proteins into the focal plane to give a more accurate representation of the 2D system (Fig. 3b, c). Further accuracy analysis can be found in Supplementary: Accuracy of the 1D-FDM and 1D-uFDM when simulating FRAP.
As FRAP is used to estimate diffusion coefficients we performed FRAP analysis on our simulated FRAP data, "Methods: Illustrative example 2: Parameters and FRAP analysis". The FRAP recovery curves (Fig. 3d) show the mean florescence recovery within the region of interest (ROI, dotted lines Fig. 3a). The results of FRAP analysis can be found in Table 1. The half time, t 1=2 seconds, is the time needed for the florescence intensity to reach half its maximum recovery. For the 2D system we calculated the half time and estimated the diffusion coefficient, e D μm 2 s À1 , using the full 2D solution and the data along a slice through the centre of the bleached area. Both methods provided comparable results, estimating the diffusion coefficient accurately to one decimal place (Table 1). FRAP analysis on the 1D-FDM estimated the diffusion coefficient to be half its actual value. In order for a 1D-FDM solution to estimate a diffusion coefficient e D ¼ 0:1 μm 2 s À1 , the actual diffusion coefficient had to be increased to D ¼ 0:22 μm 2 s À1 (Table 1). FRAP analysis on the 1D-uFDM solution resulted in an accurate estimation of the diffusion coefficient. Thus, the 1D-uFDM is able to estimate 2D diffusive movement into a trough. Furthermore, these results show that data obtained from membrane FRAP experiments could be used to fit parameters in 1D-uFDM reduced-dimension models.
Illustrative example 3: Reaction-diffusion dynamics. Spatial models rarely focus solely on diffusion. We asked to what extent could 1D-uFDM reaction-diffusion (RD) model capture 2D RD dynamics along a slide through the full system's symmetry. To address this question, a two component, mass conserved, substrate depletion model was used 8 (Fig. 4a, Methods: Illustrative example 3: RD parameters and simulations).
∇ 2 denotes diffusion, the Laplacian, it is this part of the RD equation which will be estimated by the finite difference methods ( Eqs. 1-3). Initial concentrations were set to represent a signalling event which caused a pulse conversion of molecules to u 1 from the u 2 pool (Fig. 4b). After this initial signalling event the dynamics of the system were described by the RD equations. As in previous illustrative examples, the RD equations were solved in the full 2D system and the two 1D reduced-dimension models. Diffusion was solved using the 2D-FDM, the 1D-FDM and the 1D-uFDM. The RD model was used to explore the effect of a progressively stronger initial signalling pulse on the dynamics of each solution. Results: For all solutions, larger initial signalling pulses deplete local u 2 such that the RD positive feedback becomes ineffective (u 2 1 u 2 is very small) and a u 1 trough is soon formed at the location of the initial pulse (Fig. 4c, d). Similar to the FRAP analysis results, u 1 molecules were replenished more slowly in the 1D-FDM RD solution than in the 2D-FDM and 1D-uFDM RD solutions. Thus, in the 1D-FDM RD solution, as the initial signalling pulse increases, creating a larger trough, movement into the trough becomes insufficient to reach the centre and two narrow peaks are formed (Fig. 4c, e). The 1D-FDM solution is qualitatively different from the 2D-FDM solution for α ≥ 1:5 ( Fig. 4c-e). The 1D-uFDM reproduced all the 2D-FDM RD molecular dynamics through the focal plane (Fig. 4c-e), showing that it can be used to increase the accuracy of reduced-dimension RD models. Further results and accuracy analysis can be found in Supplementary: Reaction-diffusion comparisons and the effect of geometry.

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 J of the 2D mesh, to estimate concentrations at mesh points in rows J ± 1 of the 2D mesh (Fig. 1d, f, 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 i ¼ N=2. To generate the mesh to be interpolated from re-index the 1D mesh using n ¼ i À N=2. The transformed index, to be interpolated from is To generate the mesh to be interpolated to use the transformed index to calculate the distance of each point in rows J ± 1 from the centre of the patch using the equation (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, n ¼ 0; 1; The accuracy of estimating concentrations in rows J ± 1 of the 2D mesh, using interpolation, is discussed in Supplementary: Accuracy of estimating concentrations at phantom points using interpolation.
Explicit 1D-uFDM derivation. The explicit 2D-FDM is where d x ¼ 4t 4x 2 D and d y ¼ 4t 4y 2 D. J 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 J. Thus, we have to estimate the concentrations u τ i;JÀ1 and u τ i;Jþ1 . To achieve this we use the property of radial symmetry exhibited by a patch of proteins. Transform the 4x mesh points such that the centre of the patch is at mesh point n ¼ 0 (see section "Methods: Generating the interpolation mesh in 1D"). The property of radial symmetry dictates that, for k ≥ 0, u τ k;Jþ1 ¼ u τ k;JÀ1 ¼ u τ Àk;Jþ1 ¼ u τ Àk;JÀ1 (Fig. 1d, f). Denote the estimated concentration u τ n;J ± 1 as u τ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi (Fig. 1g).
Substituting the interpolated values of u τ n;J ± 1 into the 2D-FDM equations, and removing the J subscript we get the 1D-uFDM, Explicit 1D-uFDM solution. To solve the explicit 1D-uFDM we write the 1D-uFDM equation in matrix form where u τ denotes the N × 1 vector of concentrations u τ n on mesh points ð9Þ a tridiagonal N × N matrix, with periodic boundary conditions. Using interpolation on the matrix form we derive the solution to the explicit 1D-uFDM, Explicit 1D-uFDM numerical stability condition. Let λ denote the vector of eigenvalues for matrix A. To calculate the stability condition for 1D-uFDM numerical stability recall that the valuesũ τ j are interpolated from u τ at time τ, thus, as long as a stable interpolation method is used, the solution will be stable if jλj ≤ 1 8 λ 2 λ. Gerschgorin's circle theorem 21 states that λ is bounded by the inequality, which can be simplified to, Thus, the 1D-uFDM solution will be numerically stable if 1 À 2d y ≤ 1 and À1 ≤ 1 À 4d x À 2d y . 1 À 2d y ≤ 1 is always satisfied as d y > 0. The inequality À1 ≤ 1 À 4d x À 2d y leads to the explicit 1D-uFDM stability condition, The explicit 1D-uFDM stability condition is numerically verified in Supplementary: Numerical verification of stability conditions.
Semi-implicit 1D-uFDM: derivation. A fully implicit 1D-uFDM is ill defined as molecular movement through the focal plane is inferred using the concentrations on the focal plane, Supplementary: An implicit 1D-uFDM is ill defined. A semiimplicit numerical solver can be defined in which molecular movement through the focal plane is solved explicitly and molecular movement on the focal plane is solved implicitly. The semi-implicit 2D-FDM equation is Using the same reasoning used for the derivation of the explicit 1D-uFDM, the  semi-implicit 1D-uFDM is Note, the semi-implicit 1D-uFDM scheme mirrors the 1D Crank-Nicolson method which converges and is unconditionally stable 21 . 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, a tridiagonal N × N matrix, with periodic boundary conditions. Again, using interpolation on the matrix form of the equation we derive the solution to the semi-implicit 1D-uFDM, Semi-implicit 1D-uFDM: numerical stability condition. Let λ be the vector of eigenvalues for C. For the solution of the semi-implicit 1D-uFDM to be stable two inequalities must hold, 1=λ ≤ 1 and 1 À 2d y ≤ 1. The first inequality 1 ≤ λ j j can be investigated using Gerschgorin's circle theorem 21 . Gerschgorin's circle theorem states that which can be simplified to Thus, the inequality 1 ≤ jλj is always satisfied. The second inequality j1 À 2d y j ≤ 1 expands to À1 ≤ 1 À 2d y ≤ 1. 1 À 2d y ≤ 1 as d y ≥ 0. À1 ≤ 1 À 2d y leads to the semiimplicit 1D-uFDM stability condition, The semi-implicit 1D-uFDM stability condition is numerically verified in Supplementary: Numerical verification of stability conditions. Illustrative example 1: Parameters. For results shown in Fig. 2b-e, 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: D ¼ 0:1 μm 2 s À1 , Δx ¼ Δy ¼ 0:1 μm, Δt ¼ 0:01s. These values were chosen to satisfy the numerical stability conditions of all three numerical methods. The explicit 1D-uFDM accuracy analysis found in Supplementary: Steady state accuracy of 1D-uFDMs and Accuracy dynamics of 1D-uFDMs was also taken into consideration.
Illustrative Example 2: Parameters and FRAP analysis. For the comparative FRAP analysis, D ¼ 0:1 μm 2 s À1 , Δx ¼ Δy ¼ 0:1 μm, Δt ¼ 0:01s for all solutions. The 2D ROI used to calculate the FRAP recovery curve was a circle radius 0:55 μm and the 1D ROI a line 1:1 μm in length (Fig. 3a). Both ROIs were placed in the centre of the bleached region. t 1=2 was calculated using linear interpolation on the Fig. 4 Illustrative example 3: Reaction-diffusion initial conditions and results. a Cartoon of the RD model. Black arrows represent reactions and green dashed arrows represent diffusion. Mass is conserved, molecules are only changed, never created or lost. b Initial concentrations of u 1 and u 2 for the RD models. α ¼ ½0:5; 1; ; 3 is the maximum concentration of the u 1 peak after the initial signalling pulse. u 1 x; y; 0 ð Þ¼αe À x 2 þy 2 ð Þ , u 1 x; 0 ð Þ¼αe Àx 2 and u 2 ¼ m À u 1 , where m is the mass of the system, see "Methods: Illustrative example 3: RD parameters and simulations". c Kymographs of the u 1 concentration in the 2D solution along a slice through the centre of the initial signalling pulse, and the 1D solutions. d u 1 concentration dynamics at the centre of the initial signalling pulse in the 1D and 2D solutions. e Steady-state colorplots of the 2D and 1D u 1 solutions.
FRAP recovery curve (Fig. 3d). To estimate the diffusion coefficient the equation, was used 22 , where r e represents the value of the effective radius and r n the laser radius. For the FRAP simulations, r e ¼ p 2 and r n ¼ 0:55.
Illustrative example 3: RD parameters and simulations. Symmetry breaking parameters, and a mass sufficient for the system to exhibit saturation due to substrate depletion (m ¼ 3) were chosen 8 . The mass of the system was defined as the mean concentration of molecules, where, for k ¼ 1; 2 f g, for the 2D case, and, for 1D. Δx ¼ Δy ¼ 0:1 μm. Solving the RD equations was a done using a two-step process: reactions were solved using Euler's method at time steps Δt R ¼ 1 10 À5 s for all simulations, diffusion was solved using explicit 2D-FDM, 1D-FDM and 1D-uFDM at time steps Δt D ¼ 0:002s to ensure numerical stability of the explicit 2D-FDM.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The underlying data are present in the manuscript itself and its supplementary information file. Any other remaining information can be obtained from the corresponding author upon reasonable request.